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ABSTRACT 

This  thesis  consists  of  an  interactive  program  that  enables  the  student  to  study  the 
orbital  motion  of  satellites  around  the  earth.  The  student  can  investigate  the  shape  of 
a  variety  of  orbits  by  varying  the  initial  position  and  velocity  of  the  satellite,  or  by  sup- 
plying select  orbital  parameters  i.e.  initial  orbital  radius,  eccentricity,  and  inclination. 
Satellite  maneuvers  can  also  be  studied,  like  transfer  orbits  and  inclination  changes,  by 
command  velocity  changes  at  any  location  in  the  orbit.  Also  the  effects  of  the  perturb- 
ing forces  due  to  the  oblateness  of  the  earth,  drag  for  low  earth  orbits,  and  gravitational 
attraction  from  the  sun  and  moon  can  be  investigated.  The  orbits  are  displayed  in  either 
the  perifocal  coordinate  system  around  a  model  of  the  earth,  or  the  ground  track  can 
be  displayed  on  a  map  of  the  world.  Orbital  data  is  displayed  below  the  orbital  plot. 
The  display  is  enabled  by  the  use  of  display  integrated  software  system  and  plotting 
lanauase  (DISSPLA)  subroutines. 
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I.     INTRODUCTION 

A  visual  aid  lor  students  new  to  orbital  mechanics  is  required  to  comprehend  fully 
the  dynamics  of  orbital  motion.  This  program  is  an  interactive  time  step  simulation 
program  that  calculates  and  plots  either  unperturbed  or  perturbed  elliptical  orbits.  The 
program  interacts  with  the  student  in  developing  the  initial  orbit.  Also  the  program 
enables  the  student  with  the  ability  to  change  the  velocity  of  the  satellite  at  a  specific 
location  in  the  orbit.  This  feature  will  permit  the  student  to  investigate  the  effects  of 
commanded  velocity  changes  as  in  perigee  kicks,  apogee  kicks  and  inclination  changes. 
The  user  can  also  modify  the  initial  position  and  velocity  of  the  satellite  at  the  com- 
pletion of  any  orbit. 

The  student  is  given  an  opportunity  to  investigate  the  effects  of  perturbing  forces 
on  the  satellites  orbit  by  choosing  to  have  the  program  calculate  the  orbit  with  or  with- 
out perturbing  forces.  The  variation  of  parameters  method,  as  seen  in  [Ref.  1:  pp. 
396-407],  is  used  in  calculating  the  perturbing  orbit.  The  perturbing  forces  taken  into 
consideration  are  the  following: 

1.  the  oblateness  of  the  earth 

2.  drag  for  low  earth  orbits 

3.  gravitational  force  of  the  moon 

4.  gravitational  force  of  the  sun 

In  order  to  review  fully  the  operation  of  the  program  (included  in  appendix  A)  and 
to  uncover  any  problems  or  limitations  that  plagued  the  programming,  the  program  has 
been  divided  up  as  follows: 

1.  program  design 

2.  unperturbed  orbit 

3.  perturbed  orbit 

4.  velocity  changes 

5.  graphical  plots 

The  programming  approach  and  equations  used  in  each  of  the  above  sections  will  be 
examined  in  there  respective  chapters.  A  review  of  the  coordinate  systems  used  and  their 


transformations  between  them  are  included  in  appendix  B.  Since  all  the  equations  used 
in  the  calculation  of  the  orbital  elements  are  from  reference  1.  they  will  not  be  reviewed 
in  each  chapter  but  will  be  included  in  appendix  C  for  a  quick  reference.  Equations  from 
other  sources  will  be  referenced  in  their  respective  chapters. 

Examples  of  perturbed  and  unperturbed  orbital  plots  for  a  variety  of  initial  orbital 
parameters  are  included  in  appendix  D.  Included  are  plots  of  low  earth  orbits,  transfer 
orbits  and  eeosvnchronous  orbits. 


II.     PROGRAM  DESIGN 

In  designing  this  program  an  attempt  was  made  to  make  it  not  only  as  user  friendly 
as  possible,  but  also  to  make  the  program  as  simple  as  possible  to  understand.  To 
achieve  these  goals,  the  program  would  have  to  be  written  in  a  logical  manner,  in  a 
computer  language  that  is  easy  to  follow,  the  program  would  have  to  run  on  terminals 
readily  available  to  students  (at  the  Naval  Postgraduate  School  (NFS)),  and  the  program 
would  have  to  be  easily  used  by  students  with  a  minimum  amount  of  computer  or  orbital 
mechanics  knowledge. 

FORTRAN  was  chosen  as  the  programming  language  since  it  is  a  wildly  used  sci- 
entific language  and  it  allows  for  very  structured  programming.  By  programming  in  a 
structured  format,  the  program  can  be  expanded  in  the  future  with  a  minimum  amount 
of  time  required  to  understand  the  programming  code.  FORTRAN"  also  allows  for 
double  precision  numbers  to  be  used  in  the  calculation  of  the  orbit.  This  is  critical  when 
round  off  error  in  single  precision  could  be  greater  then  the  actual  change  that  one  is 
trying  to  model.  The  equations  in  the  descriptions  of  the  program  might  not  exactly 
match  the  equations  in  the  listings  because  of  special  programming  techniques  which 
must  be  included  in  most  computer  programs  to  handle  such  problems  as  "division  by 
zero ". 

The  display  integrated  software  system  and  plotting  language  (DISSPLA)  package 
available  on  the  mainframe  computer  at  NTS  was  used  to  enable  a  variety  of  graphical 
displays  with  a  minimum  amount  of  programming.  DISSPLA  has  a  set  of  subroutines 
that  the  programmer  calls  to  display  data  contained  in  arrays.  This  requirement  forces 
the  program  to  load  arrays  with  the  satellites  position  in  order  for  it  to  be  plotted.  The 
TEC618  computer  terminal  and  associative  plotter  was  used  for  ease  of  gaining  hard 
copy  plots  of  the  orbits  and  the  diversity  of  locations  that  are  available  here  at  NTS. 
In  order  to  run  a  program  in  DISSPLA  the  user  must  first  define  storage  space  of  1500k 
and  designate  temporary  disk  space,  and  then  call  DISSPLA  with  the  program  name. 
This  is  accomplished  with  the  following  commands: 
1.   DEFINE  STORAGE  1500K 


:.  i  cms 

3.  TDISK  4  DIS 

4.  DISSPLA  ORBIT 

To  make  the  program  user  friendly,  the  user  is  prompted  for  inputs  via  the  keyboard. 
The  entry  is  usually  a  number.  A  yes  or  no  response  can  be  entered  by  typing  "Y"  or  a 
"X".  In  most  cases  the  program  does  a  check  to  see  if  the  input  is  appropriate.  In  order 
to  make  it  as  easy  as  possible  for  the  student  to  get  the  desired  orbit  displayed,  the 
program  requires  only  the  initial  position  and  velocity  of  the  satellite.  The  initial  posi- 
tion and  velocity  of  the  satellite  is  supplied  by  the  user  in  one  of  two  ways.  The  user  can 
input  the  position  and  velocity  of  the  satellite,  using  the  perifocal  coordinate  system 
(UK),  or  the  user  can  let  the  program  place  the  satellite  on  the  "I"  axis  of  the  UK  system 
at  the  radius  of  perigee  (RP)  distance  supplied  by  the  user.  This  latter  choice  gives  the 
initial  location  of  the  satellite,  but  to  get  the  velocity  the  program  will  prompt  the  user 
for  one  of  the  following: 

1.  the  actual  velocity  in  the  UK  system. 

2.  the  eccentricity  (e)  of  the  orbit.  In  which  case  the  velocity  is  calculated  from  the 
following  equations: 

RP 

a  = =  semi-major  axis 

1  —  c 

u 
EMl  =  —  -^—  =  energy  mass 

Where  a  =  MG 

M  =  mass  of  earth 

G  =  Universal  eravitational  constant 


u 


3.   the  radius  of  apogee  (RA).  The  velocity  is  calculated  by  first  calculating  the  ec- 
centricity (e)  from  the  following: 

RA  -  RP 

e  = 


RA  -r  RP 


With  the  eccentricity  the  same  equations  used  above  are  used  to  calculate  the  ve- 
locity. 

In  order  to  give  the  velocity  a  direction  the  inclination  (i)  of  the  orbit  is  required  from 
the  user.   The  following  equations  are  used  to  calculate  the  velocity  vector: 

r,  =  0.0 


Vj-  vcos(i) 
vk  =  v  sin(/) 

The  program  will  check  to  ensure  that  the  orbital  eccentricity  is  less  than  1.0,  if  it  is  not 
then  the  program  will  reject  the  inputs.  After  the  initial  input  are  accepted,  the  program 
will  do  calculations  for  the  six  orbital  elements  required  to  describe  the  size,  shape  and 
orientation  of  the  orbit,  and  to  pinpoint  the  position  of  the  satellite  along  the  orbit  at  a 
particular  time.    This  classical  set  of  six  orbital  elements  are  as  follows: 

1.  a,  semi-major  axis. 

2.  e.  eccentricity. 

3.  i,  inclination. 

4.  Q.  longitude  of  the  ascending  node. 

5.  co.  argument  of  perigee  passage. 

6.  T. time  of  perigee  passage. 

The  program  actually  calculates  more  orbital  elements  than  the  six  classical  elements 
required  to  plot  the  orbit,  this  is  done  in  an  effort  to  make  the  program  as  robust  as 
possible.   This  will  add  in  the  ability  to  expand  the  program  in  the  future. 

If  the  satellite  is  not  initially  at  the  perigee  point  then  the  satellite  must  first  be 
stepped  around  to  the  perigee  point.  The  program  then  enters  a  loop  that  calculates  the 
orbit  from  the  perigee  point  through  one  complete  orbit  around  the  earth  and  back  to 
the  perigee  point.  The  orbit  is  calculated  in  steps  of  2  times  pi  divided  by  an  integer,  i.e., 
2  times  pi  divided  by  50.  This  step  size  was  used  to  ensure  a  smooth  orbit  for  display 
purposes  and  also  to  get  within  adequate  distance  to  the  perigee  point  or  other  location 
for  a  velocity  change.  After  the  loop  is  completed,  the  program  will  offer  the  user  a 
choice  of  the  following  plots  to  check  the  orbit: 

1.  perifocal 

2.  groundtrack 

The  program  then  goes  into  a  loop  offering  the  user  the  following  choices  until  the  user 
decides  to  end  the  program: 

1.     plot  another  view  of  the  same  orbit. 
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III.     UNPERTURBED  ORBIT 

The  subroutines  that  calculate  the  unperturbed  orbit  are  the  most  widely  used  sub- 
routines in  the  entire  program.  These  subroutines  are  called  to  step  the  satellite  around 
to  the  perigee  point  from  the  user  supplied  initial  position  and  velocity,  to  calculate  the 
next  unperturbed  orbit,  and  for  any  velocity  change.  No  matter  which  of  these  sources 
supply  the  initial  position  and  velocity  the  program  calculates  the  unperturbed  orbit  in 
the  same  manner.  The  only  difference  is  where  in  the  orbit  the  satellite  is  initially  when 
these  subroutines  are  called.  Before  the  unperturbed  subroutines  are  called,  the  orbital 
elements  are  calculated. 

The  unperturbed  subroutines  are  called  by  a  single  subroutine  UNPRET'  which  has 
the  following  basic  algorithm: 

1.  Increment  time  by  the  time  step  size  (DT).  The  time  step  was  chosen  as  the  period 
divided  by  fifty  to  give  a  smooth  plot,  but  more  importantly  to  ensure  that  the 
satellite  is  within  an  acceptable  distance  from  a  specific  location  for  a  velocity 
change.  The  angular  error  caused  by  the  step  size  can  be  as  much  as  PI  50  from 
the  desired  point  for  a  circular  orbit  and  will  increase  for  more  eccentric  orbits. 
This  error  becomes  a  factor  when  the  user  is  making  velocity  changes,  and  therefore 
it  will  be  covered  in  that  chapter  in  further  detail. 

2.  Calculate  the  new  elements.  The  calculation  of  the  new  elements  is  the  heart  of  this 
algorithm.  The  size,  shape  and  orientation  of  the  orbit  remains  unchanged.  What 
is  required  is  the  position  of  the  satellite  along  the  orbit  as  a  function  of  time.  The 
problem  becomes  a  matter  to  solve  "the  Kepler  problem'-predicting  the  future  po- 
sition and  velocity  of  an  orbiting  object  as  a  function  of  some  known  initial  posi- 
tion and  velocity  and  the  time  of  flight  [Ref.  1:  p.  181].  An  algorithm  using  these 
principles  will  follow: 

a.  A  time  step  (DT)  is  added  to  the  time  of  flight(TF).  time  of  flight  is  the  elapsed 
time  since  the  satellite  passed  the  perigee  point. 

TF=  TF+DT 

b.  The  new  mean  anomaly  (MA)  is  calculated  from  the  new  time  of  flight,  and  the 
mean  motion  (MM). 

MA  =  MM  x  TF 

c.  With  the  new  mean  anomaly  the  new  eccentric  anomaly  (EA)  is  calculated. 
Because  the  solution  to  the  Kepler  problem  (MA  =  EA  —  e  x  sin(EA))  is 
transcendental,  an  iterative  solution  based  on  the  Newton  method  of  root  find- 
ing is  used.  The  root  in  question  is  a  solution  to  the  equation 
[MA  -  EA  +  e  x  s'm(EA)  =  0)  .  This  algorithm  takes  the  form  of  [Ref.  I:  p.  222]: 

1 1  MAn  =  EAn  -ex  sin(EAr) 


2) 

MA -MA, 


EAnM  =EAn  +  -  , 

'•  '  {I  -ex  cos{EA„)) 

Where  this  equation  is  applied  initially  to  EA0  =  MA  and  then  reapplied 
until  the  difference  between  MA  and  MA„  becomes  small  enough  to  be  ig- 
nored. 

d.  The  new  true  anomaly  (ve)  is  calculated  from: 

cos"  (e  —  cos(  EA)) 


vo  = 


e  cos{EA)  —  1 


3.   Calculate  the  new  position  and  velocity.    The  position  and  velocity  are  calculated 
in  the  perifocal  coordinate  system  (PQW).    The  PQW  system  uses  the  orbit  as  its 

fundamental  plane  and  therefore  requires  only  two  coordinate  to  specify  the  satel- 
lite's position  and  velocity.  The  zrt  coordinate  is  by  definition  always  equal  to  zero. 
The  position  of  the  satellite  is  calculated  as: 


xw  =  r  cos  v 

yw  =  r  sin  v 

-,  =  o 

The  velocitv  of  the  satellite  is  calculated  as: 


-,.  =  v  jr  (  -  sin  v0) 

'  u 
v  =     /—  (e  +  cosv0) 


\    P 


A.  Store  position  and  elements  in  arrays  for  plotting.  In  order  for  the  program  to  plot 
the  orbit  the  radius,  true  anomaly,  inclination,  and  argument  of  perigee  must  be 
stored  in  arrays.  The  use  of  these  arrays  to  plot  the  orbit  will  be  explained  in 
chapter  6. 

5.  The  process  is  repeated  until  the  satellite  is  at  the  perigee  point  and  the  true 
anomaly  is  two  pi. 

The  procedure  used  to  calculate  the  unperturbed  orbit  leave  very  little  to  be  modified 
by  a  programmer.  The  only  choices  that  had  to  be  made  concerned  step  size,  how  to  tell 
the  UNPRET  subroutine  that  the  perigee  point  had  been  reached,  and  a  value  of  ac- 
ceptable error  for  newtons  method.  For  the  unperturbed  orbit,  the  step  size  just  had  to 
be  small  enough  to  produce  a  smooth  plot  of  the  orbit.  Two  indicators  for  perigee  were 
used,  one  was  that  the  true  anomaly  was  greater  than  6.21  radians  (two  pi  equals  6.2S 
radians)  and  the  time  from  the  previous  perigee  point  will  be  greater  then  the  period. 
The  two  indicators  were  logically  'and'  together  to  ensure  the  perigee  point  was  reached. 


The  disparity  between  two  pi  and  0.21  radians  is  due  to  the  error  produced  by  the  sat- 
ellite not  beginning  the  orbit  at  exactly  the  perigee  point  and  the  step  size  used  go 
around  the  orbit.  The  acceptable  size  of  error  for  newtons  method  was  set  at  l.OE-10, 
because  for  an  unperturbed  orbit  this  would  be  the  major  contributor  to  any  error  in  the 
orbit  and  the  magnitude  of  this  error  would  be  acceptable.  However:  in  a  perturbed 
orbit  there  are  other  factors  contributing  to  determining  the  acceptable  error,  and  these 
will  be  discussed  in  the  next  chapter. 


IV.     PERTURBED  ORBIT 

The  perturbed  orbit  uses  the  same  basic  routines  as  the  unperturbed  orbit  in  step- 
ping the  satellite  around  the  earth  with  one  major  difference,  the  perturbing  forces 
produce  a  time  rate  of  change  of  the  orbital  elements  that  must  be  applied  at  each  time 
step.  The  variation  of  parameters  method  is  used  to  determine  this  influence  of  the 
perturbing  forces  on  the  orbital  elements.  The  analysis  is  simplified  by  using  the  orbital 
coordinate  system  'RSW,  as  explained  in  appendix  B.  The  basic  algorithm  is  as  follows 
[Ref.  I:  p.  407]: 

1.  At  i  =  t0  calculate  six  orbital  elements. 

2.  Compute  the  perturbing  forces  and  transform  it  at  i  —  i0  to  the   RSW  SYSTEM. 

3.  Compute  the  time  rate-of-change  of  the  elements. 

-4.   Calculate  the  change  of  elements  for  one  time  step,  and  add  the  changes  to  the  old 

values  at  each  step  to  get  the  new  elements. 

5.  From  the  new  values  of  the  orbital  elements,  calculate  a  position  and  velocity. 

6.  Go  to  the  step  2  and  repeat  until  the  final  time  is  reached. 

The  steps  in  the  algorithm  will  be  explained  in  the  following  sections: 

A.     ORBITAL  ELEMENTS 

The  standard  orbital  elements  a.  c,  i.  Q,  co  and  T  (or  M)  will  be  used,  where 

a  =  semi-major  axis 

e  =  eccentricity 

i  =  inclination 

Q  =  longitude  of  ascending  node 

a)  =  argument  of  perigee 

T  =  time  of  perigee  passage 

(\I0  =  mean  anomaly  at  epoch  =  A/  —  n(t  —  ta)).  The  elements  are  calculated  only 
at  the  beginning  of  the  orbit  from  the  initial  position  and  velocity  vectors.  The  elements 
are  then  changed  continuously  throughout  the  orbit  by  adding  the  changes  due  to  the 
perturbing  forces.  For  the  perturbed  orbit,  the  satellite  will  always  begin  at  the  perigee 
point.    This  is  done  so  one  complete  orbit  is  from  perigee  point  to  perigee  point. 


B.     COMPUTE  PERTURBING  EORCES 

The  variation  of  parameters  method  requires  that  the  perturbing  forces  be  calculated  at 
each  step  in  the  orbit.  In  order  to  do  this  a  model  of  each  perturbing  force  must  be 
developed.  The  following  perturbing  forces  where  used  in  calculating  the  total  perturb- 
ing force  effecting  the  satellite: 

1.  oblateness  of  the  earth 

2.  atmospheric  drag 

3.  gravitational  attraction  of  the  sun 

4.  gravitational  attraction  of  the  moon 

The  magnitudes  of  these  forces  have  an  enormous  range  of  values  and  are  dependent 
on  the  distance  the  satellite  is  from  the  perturbing  body.  Figure  1  on  page  12  shows  a 
graphical  representation  of  the  magnitude  of  the  perturbing  forces  in  a  log-log  plot  of 
perturbing  forces  per  unit  mass  [Ref.  2:  p.  IV-61].  The  model  of  each  of  these  forces 
follows: 

1.     NON-SPHERICAL  EARTH 

The  earth  is  not  perfectly  spherical,  but  bulges  around  the  equator.  The  polar 
and  equatorial  diameters  are  12713.0  Km  and  12756.3  Km.  respectively.  The  oblateness 
results  in  a  perturbing  force  per  unit  mass  with  these  components  in  the  'RSW  coordi- 
nate system  [Ref.  3:  p.  SI]: 

I  —}uJ-r~)  2  -, 

Fr  = f —  (1  —  3  sin  (/)  sin"(w0)) 


Fs  = : (  sin  (/)  sin(r/0)  cos(>0)) 


{  —  3u7-,r!) 

Fw  = r^^  ( sm(/)  cos(/)  sin(w0)) 

r 

The  variable  and  constants  of  these  equations  are  defined  below: 

1.   Variables: 

a.   u0  =  the  argument  of  latitude  and  is  equal  to  the  true  anomaly  v0  plus  the  ar- 
gument of  perigee  cj. 


b.   r  =  the  radius  from  the  center  of  the  earth  to  the  satellite. 
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Figure  1. 


Comparison  of  perturbation  magnitudes. 
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r=\r  i 
2.   Constants: 

a.   u  =  the  gravitational  parameter  of  the  earth, 


Ai  =  398601.2 


(km2) 


b.  J2  =  the  second  harmonic  of  oblateness  coefficient,  determined  by  experimental 
observations, 

J2=  1.0823£-3 

c.  rt  =  the  mean  radius  of  the  earth, 
re  =  631S2E3Km 

2.     ATMOSPHERIC  DRAG 

The  formulation  of  atmospheric  drag  equations  are  plagued  with  uncertainties 
of  atmospheric  fluctuations,  frontal  areas  of  orbiting  object  (if  not  constant),  the  drag 
coefficient,  and  other  parameters.  A  fairly  simple  formulation  will  be  given  here.  Drag, 
by  definition,  will  be  opposite  to  the  velocity  of  the  vehicle  relative  to  the  atmosphere. 
Thus,  the  perturbing  force  is 

F  =  -{  -T-  ).CD*AR.  DEX  .  v .  v 

2m 

The  velocity  vector  is  in  the  'UK'  system  so  the  resulting  force  is  also  in  the    UK'  sys- 
tem.  Therefore  a  transformation  to  the  'RSW  system  is  required. 

The  variables  and  constants  of  this  equation  are  defined  below: 
1.   Variables: 

a.  v  =  speed  of  vehicle. 

b.  CD  =  the  dimensionless  drag  coefficient.  The  drag  coefficient  CD  has  a  value 
between  1  and  2.  It  takes  a  value  near  1  when  the  mean  free  path  of  the  at- 
mospheric molecules  is  small  compared  with  the  satellite  size,  and  takes  a  value 
close  to  2  when  the  mean  free  path  is  large  compared  with  the  size  of  the  satel- 
lite. The  drag  coefficient  will  be  modeled  with  CD  =  2  when  the  satellites  alti- 
tude is  greater  than  550km  and  equal  to  1  otherwise.  [Ref  4:  p.  295] 

c.  DEN  =  atmospheric  density  at  the  vehicle's  altitude.  The  density  is  spherically 
symmetric,  and  will  be  modeled  using  exponential  steps  using  the  parameters  in 
Table   1  on  page  14  and  the  following  formula  [Ref.  1:  pp.  423-424]: 

3{z)  =^v<-'z>  . 
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Table  1 

.     ATMOSPHERIC  PARAMETERS  AND  ^ 

ALl  ES 

/■  km) 

|  k 

j 

5{z) 

0-150 

1.225E-02 

4.74E-02              !  en 

1  ?">2^E-02 

150 

1.0E-03 

150-550 

1.79846E-01 

4.3614E-02 

550 

3.0E-8 

550 
> 

1.015484E-07 

2.2169SE-07 

1 500 

3.65E-09 

4100 

1.0E-12 

2.   Constants  set  to  typical  values: 

a.  m  =  mass  of  the  satellite,  set  equal  to  100kg. 

b.  AR  =   the  cross-sectional  area  of  the  vehicle  perpendicular  to  the  direction  of 
motion,  set  equal  to  20m: 

3.  PERTURBING  FORCE  DUE  TO  HEAVENLY  BODY 

The  satellite  will  experience  perturbation  forces  due  to  the  gravitational  effects 
of  the  sun  and  the  moon.  The  perturbation  force  from  a  perturbing  body  is  the  differ- 
ence between  the  gravitational  force  due  to  the  perturbing  body  at  the  satellite  and  the 
gravitational  force  the  satellite  would  experience  if  it  were  at  the  center  of  the  earth. 
From  Figure  2  on  page  15.  the  perturbing  force  per  unit  mass  of  the  satellite  is 


U-r, 


rnir,  —  rir 


':JP  ~  >'lr 


/V 


F  P 


r, 


The  variable  and  constants  are  defined  below: 

1.  Variables: 

a.  rp  =  distance  from  the  earth  center  for  the  perturbing  body 

b.  ip  =  unit  vector  from  the  earth  to  the  perturbing  body 

c.  r    =  distance  from  earth  center  to  the  satellite 

d.  i,  =  unit  vector  from  the  earth  to  the  satellite 

2.  Constants: 

a.   \jl.  =  gravitational  constant  of  the  perturbing  body  =  MFG 

The  subscript  p  is  to  be  replaced  by  s  if  the  perturbing  body  is  the  sun.  and  by  m  if  the 
perturbing  body  is  the  moon.  We  will  assume  that  r  <  <  rF  then  the  equation  above 
bvxo:..o^ 
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Figure  2. 


Perturbation  forces. 


Fp  =  ( — )( 7~  )(3(i>/p)ip  -  /r) 


The  unit  vectors  /,  and  /,  can  be  written  in  terms  of  the  'UK'  system  as: 
ir  =  (  cos(Q)  cos(uo)  -  sin(Q)  cos(/)  sin(«o))/  +  (  cos(z^)  sin(Q)  +  cos(a>)  cos(/)  sin(^))J  + 

(  sin(/)  sin(^))A' 

i  =  (  cos(Q  )  cos(uo  )  -  sin(Q-)  cos{L)  sin(uoP))7  + 


(  cos(uop)  sin(Q  )  +  cos(co p)  cos(/_)  sin(w0  ))J  +  ( sin(/,)  sin(w0  ))K 


where  Q.,  i,  and  Up  are  the  orbital  elements  of  the  satellites  and  Qn,  ip,  and  u0r  are  the 
orbital  elements  of  the  perturbing  body.  The  formulas  above  use  the  'UK'  system,  and 
as  such  the  resultant  forces  must  be  transformed  to  the  'RSW  system.  Models  of  the 
sun  and  moon  orbits  are  required  to  calculate  ~rp  and  ip.  The  models  used  in  the  program 
for  the  sun  and  moon's  orbits  follows:  [Ref.  3:  pp.  73-74] 
a.     SUN'S  POSITION 

In  order  to  model  the  suns  orbit,  a  number  of  simplifications  had  to  be 
made  in  the  actual  parameters  of  the  suns  orbit.  First  the  sun  will  be  assumed  to  be  in 
a  circular  orbit. This  means  that  the  radius  (r)  to  the  sun  will  be  constant,  and  the  ec- 
centricity (e)  will  equal  0.0  instead  of  its  true  value  of  0.017.    The  other  assumption  will 
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be  to  place  the  sun  on  the  T  axis  of  the  UK'  system  at  the  beginning  of  the  program 
and  have  it  progress  through  its  orbit  as  the  program  runs.  These  changes  will  not  effect 
the  perturbing  force  in  any  noticeable  magnitude. 

The  following  variables  and  constants  where  used  in  the  program  to  model 
the  suns  orbit  after  applying  the  simplifications:    [Ref.  3:  pp.  75-78] 
1.   Constants: 

(Am2) 

a.  Gravitational  Constant:  0  =  6.6/ £—  11  — — — 

kg2 

b.  Sun's  Mass:  m  =  1.99£30A> 


c.    Sun's  Gravitational  parameter: 
/i,=  1.32733£20Jy 


Am 


d.  Sun's  eccentricity:  <?.  =  0.0 

e.  Radius  of  orbit,  assume  sun  is  in  circular  orbit:  rs  =  1.49£ll/w 

f.  Sun's  inclination:    si  =  23.45  deg.  =  4.09279709d-01  radians 

g.  Longitude  of  ascending  node:  Q,  =  0.0 
h.  Argument  of  perigee:  cot  =  0.0 

2.   Variables: 

a.  The  true  anomaly  of  the  sun's  position  as  a  function  of  the  time  the  satellite  has 
been  in  orbit: 

->_ 

'    AfT)=  356  x  2*4  x  3600"  ^ 
Where  TT  =  true  time,  the  time  the  satellite  has  been  in  orbit  (sec) 

b.  Sun's  Position  vector:  r  —  r  cos  v0sP  +  r  sin  v0;O 

r, 

c.  Lnit  vector  from  the  earth  to  the  sun:  i.  =  — r— 

'        I  r,  | 

b.     MOON'S  POSITION 

In  modeling  the  orbit  of  the  moon,  similar  assumptions  where  used  as  with 
the  sun.  The  moons  orbit  will  be  assumed  to  be  circular,  actually  the  eccentricity  is 
equal  to  0.055.  By  placing  the  moon  initially  on  the  T  axis  of  the  'UK'  system  along 
with  the  sun.  the  gravitational  forces  of  the  two  bodies  will  combine  to  a  maximum. 
However;  since  the  moons  orbital  period  is  only  27.3  days,  the  moon  will  not  stay  in  this 
alignment  and  the  magnitude  of  the  combined  forces  will  van'  with  time.  The  inclination 
of  the  moons  orbit  is  not  constant,  but  drifts  between  IS. 3  and  2S.6  decrees  in  ten  vears. 
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Also  the  longitude  of  the  ascending  node  (CI)  oscillates  between  13  and  -13  degrees.  To 
simplify  this  the  inclination  will  be  chosen  as  a  constant  23.5  degrees  and  the  longitude 
of  the  ascending  node  as  0.0  degrees.  For  the  time  period  involved  in  calculating  the 
perturbed  orbit,  these  assumptions  will  not  make  any  significant  difference. 

The  following  variables  and  constants  were  used  in  the  program  to  model 
the  moons  orbit,  after  applying  the  simplifications: 
1.   Constants: 

(Xm2) 

a.  Gravitational  Constant:  G  =  6.67 E  —  1 1  — — — 

kg2 

b.  Moon's  Mass:  mm  =  7.35£22/cg 

(Xm2) 

c.  Moon  s  Gravitational  Parameter:  \xm  —  OMm  =  4.90E12 


d.  Moon's  eccentricity:  em  =  0.0 

e.  Radius  of  orbit,  assume  moon  is  in  circular  orbit:  rm  =  3.844£8&m 

f.  Moon's  inclination:  i  =  23.5deg.  =  4.10152374E-1  radians 

g.  Moon's  longitude  of  ascending  node:  Qm  =  0.0 

h.  Moon's  argument  of  perigee:  u„  =  0.0 

i.  Moon's  period:  T=  27.3  days  [period] 
2.   Variables: 

a.  The  true  anomaly  of  the  moon's  position  as  a  function  of  the  time  the  satellite 

->_ 

has  been  in  orbit:  \\JTT)  =-rzri — rr — — —  7T 

2  .3  x  24  x  j 600 

b.  Moon's  position  Vector:  r  =  r  cos  vBmP  +  r  sin  v0mQ 

c.  L  nit  vector  from  earth  to  moon:  im  =  — - — 

I    I'm    I 

The  models  of  the  sun  and  moons  orbit  calculates  the  position  vector  in  the  TQW 
system  and  therefore  the  position  vector  must  be  transformed  to  the  'UK'  system. 

C.     RATE-OF-CHANGE  OF  ORBITAL  ELEMENTS 

The  derivations  and  equations  of  the  rates-of-change  of  the  orbital  elements  are 
contained  in  reference  1  pages  398  to  406.  Therefore;  only  a  summary  of  the  actual  an- 
alytic expressions  for  the  rate-of-change  of  the  parameters  in  terms  of  the  perturbations 

will  follow: 

1.    Rate-of-change  of  the  semi-major  axis: 
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da        r      ~e  sin  vo      -,  r-    .    r   ~"\ 


/?  ^  1  —  <_• 
Where  ri  is  the  mean  motion  of  the  satellites  orbit. 


\    a 

2.   Rate-of-chanee  of  the  eccentricity; 


?.-        2, 


de  \  X  ~  e~  sinvo  \  l~e  a\\-e  , 

di  n  a  fl'a-e 

3.   Rate-of-chanse  of  the  inclination: 


at 


r  cos  z/-, 


■T— C =]F,, 

Kt  ,    2       ,  2 

k  a  -v    1  —  e 

4.  Rate-of-change  of  the  longitude  of  the  ascending  node: 
,yo  »•  sin  wn 

i7r=[ = — ]^. 

/;  a  ^  1  —  e     sin  / 

5.  Rate-of-change  of  the  argument  of  perigee: 
-  =  (^— -  )r  +  (//— -  )5  4-  (</  —  >,v. 


</;  <//  c/r  ;  <// 

Where 


do  ~  \  {  ~  e     cos  vo 

(^).-[-^][»"'-,(i+,.t.'cos,.o)if, 

,/,,,  —  >•  cot  /  sin  «n 

/?  «  %  / 1  —  e 

6.  Rate-of-change  of  the  eccentric  anomaly: 

[(  sinv0  +  --f-)(l  +<?cos  v0)  -  (  cos  v0  +  e)( -~  cos  v0  +  «  sin  v0)] 

a  LA 1 at at 

dt      '  sin(EA)  [l  +  ecosv0]2 

7.  Rate-of-change  of  the  mean  anomaly: 

dXIA        dEA        de     .  ,_,N  (E4)dEA        dn\ 

sin(/..l)  —  ex  cos ; — —  U  —  /,,) 


dt  dl  dt  dt  dt 


IS 


This    equation    reduces    to    the    following    for    circular    and    ecliptic    orbits 
(0<  =e<  1). 

dUA        -1    r  2r       (I  ~e2)  \  - /-  r  .  dn' 

— ; —  =  —  [  — -, cos  vQ]Fr  -  [  — ; ][1  H —  J  sin  v0F,  -  t  — t- 

di  ria       a  <=  uj   r      l    n  ae  a(\-e)  iU 

Where  the  Rate-ol-change  of  the  mean  motion: 

dn'  _  r    ~-\u       da 
di    ~  L  2n'a4      di 

[ref.  1  p.  396-407] 

D.     NEW  ORBITAL  ELEMENTS 

The  change  of  each  element  is  calculated  by  multiplying  the  rate-of-change  o{  the 
element  by  the  time  step  (DT).  The  change  in  the  orbital  elements  are  then  added  to  the 
current  values  of  the  elements  to  give  the  new  orbital  elements.  With  the  new  elements 
calculated,  the  satellite  is  stepped  forward  and  the  new  position  and  velocity  are  calcu- 
lated in  the  same  manner  as  the  unperturbed  orbit  (chapter  3).  Also  as  with  the  unper- 
turbed orbit,  the  process  is  repeated  until  the  satellite  is  at  the  perigee  point,  indicated 
by  the  time  of  flight  (TF)  equal  to  the  period  of  the  perturbed  orbit. 
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V.     VELOCITY  CHANGES 

The  ability  of  the  student  to  change  the  velocity  of  the  satellite  at  any  position  in 
the  orbit  is  a  vital  element  in  this  program.  With  velocity  changes  the  student  can  in- 
vestigate the  effects  of  varying  the  satellites  velocity  as  in  transfer  orbits  and  inclination 
changes.  In  order  to  simplify  the  program  the  unperturbed  orbit  is  used  throughout  this 
routine.    The  velocity  change  algorithm  used  in  the  program  follows: 

1.  Rotate  to  velocity  change  location. 

The  user  is  given  the  choice  of  changing  the  velocity  of  the  satellite  at  the 
perigee,  apogee  or  at  any  true  anomaly.  If  the  user  chooses  perigee  or  apogee  as 
the  change  locations,  the  true  anomaly  is  set  equal  to  zero  or  pi  radians  respect- 
fully. With  the  location  of  the  velocity  change,  the  satellite  is  first  stepped  around 
to  the  desired  true  anomaly.  The  stepping  is  identical  with  the  unperturbed  orbit 
with  the  exception  that  the  stepping  terminates  when  the  true  anomaly  is  greater 
or  equal  to  the  desired  true  anomaly.  With  a  step  size  of  one  fiftieth  of  the  period, 
the  satellite  is  actually  stepped  around  to  a  location  near  the  desired  location.  This 
variance  can  be  reduced  by  decreasing  the  step  size  but  this  would  increase  the 
computation  time.  This  error  will  be  a  major  factor  in  precise  calculations  of 
transfer  orbits,  or  any  other  orbital  maneuver  where  precise  velocity  changes  are 
required.  However:  this  program  is  not  a  tool  to  calculate  precise  orbital  maneu- 
vers, but  rather  a  learning  tool  for  the  student  to  get  a  feel  for  the  results  of  velocity 
changes  in  a  satellite's  orbit. 

2.  Change  the  velocity. 

With  the  satellite  at  the  desired  location,  the  program  calculates  and  displays 
for  the  user  the  satellite's  current  velocity,  escape  velocity  and  circular  velocity  (the 
velocity  required  to  circularize  the  orbit).  The  program  will  not  allow  velocities 
greater  than  or  equal  to  the  escape  velocity.  The  user  is  given  the  option  to  enter 
a  new  velocity  in  the  'UK'  system  or  to  change  the  magnitude  of  the  velocity  in  the 
orbital  plane.  If  the  user  chooses  to  change  the  velocity  in  the  orbital  plane,  the 
program  will  prompt  the  user  for  the  magnitude  of  the  velocity  change,  and  multi- 
ply this  change  by  a  unit  vector  in  the  direction  of  the  satellites' velocity.  This  ve- 
locity change  vector  is  then  added  to  the  satellites  velocity  vector,  to  calculate  the 
new  velocity  vector. 

3.  Calculate  new  elements. 

The  orbital  elements  are  calculated  with  the  new  velocity  vector  and  the  satel- 
lite's position  vector. 

4.  Complete  the  orbit. 

The  program  will  complete  the  orbit  to  the  new  perigee  point  using  the  satel- 
lite's position,  new  velocity  and  new  elements.  There  are  a  number  of  problems 
that  arise  if  the  satellite  is  just  stepped  around  to  the  perigee  point.  For  example, 
with  velocity  changes  in  the  orbital  plane  the  apogee  and  perigee  directions  can 
physically  swap.  This  is  a  problem  when  plotting  with  the  perifocal  coordinate 
system  because  the  .V,  axis  points  toward  perigee.  To  avoid  problems  like  this  the 
arrays  used  in  plotting  the  orbit  must  be  cleared  and  the  satellite's  current  position 
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and  velocity  be  treated  as  initial  conditions.  However;  to  compare  the  old  and  new 
orbits  there  is  a  desire  to  retain  as  much  of  the  previous  orbit  as  possible.  The 
velocity  changes  where  divided  into  the  following  four  cases  to  handle  these  prob- 
lems: 

a.  Change  velocity  in  the  orbital  plane  at  the  perigee  point  with  the  new  velocity 
greater  than  the  circular  velocity.  The  perigee  point  will  remain  the  same  so  the 
satellite  is  stepped  around  using  the  unperturbed  subroutines. 

b.  Change  velocity  in  the  orbital  plane  at  the  perigee  point  with  the  new  velocity 
less  than  or  equal  to  the  circular  velocity.  The  perigee  and  apogee  directions 
will  switch  so  the  plotting  arrays  are  first  cleared  and  stored  with  the  current 
location  data.  Because  the  satellite  is  now  at  the  apogee  point  the  satellite  is 
stepped  around  to  the  perigee  point  storing  the  second  half  of  the  orbit.  The 
entire  next  orbit  is  calculated  and  stored  to  get  a  complete  orbit. 

c.  Change  velocity  in  the  orbital  plane  at  the  apogee  point  with  the  new  velocity 
less  than  the  circular  velocity.  The  perigee  and  apogee  directions  will  remain  the 
same,  so  the  satellite  is  stepped  around  to  the  perigee  point  completing  the  orbit. 

d.  This  last  case  catches  all  the  following  velocity  changes;  velocity  change  in  the 
orbital  plane  at  the  apogee  point  with  the  new  velocity  greater  than  or  equal  to 
the  circular  velocity,  velocity  changes  at  any  other  true  anomaly  in  the  orbital 

plane,  and  any  velocity  change  out  of  the  orbital  plane.  The  plotting  arrays  are 
cleared  and  stored  with  the  current  location  data.  No  matter  where  in  the  orbit 
the  satellite  is.  the  satellite  is  first  stepped  around  to  the  perigee  point,  and  to 
ensure  a  complete  orbit  is  plotted  the  entire  next  orbit  is  also  calculated  and 
stored. 
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VI.     GRAPHICAL  PLOTS 

The  program  provides  two  types  of  graphical  displays  of  the  orbit,  a  display  in  the 
perifocal  coordinate  system  and  a  display  of  the  satellite's  ground  track.  Each  display 
type  is  useful  in  observing  different  aspects  of  the  orbit.  The  perifocal  display  will  allow 
the  user  to  see  how  certain  orbital  parameters  change  with  different  initial  positions  and 
velocities,  and  also  how  the  parameters  change  with  velocity  changes  at  varying  posi- 
tions in  the  orbit.  The  ground  track  will  enable  the  user  to  gain  an  appreciation  for  the 
physical  location  of  the  satellite  above  the  earth,  and  see  how  the  orbital  parameter  af- 
fects the  path  of  the  satellite.  The  ground  track  will  also  display  the  precession  of  a  se- 
quence of  orbits.  Both  displays  plot  the  position  steps  to  give  the  user  an  understanding 
of  how  the  satellite  speeds  up  at  perigee  and  slows  down  around  apogee. 

The  DISSPLA  package  on  the  mainframe  computer  was  used  to  enable  the  plotting 
of  the  orbits.  The  versatility  of  plotting  subroutines  of  DISSPLA  makes  the  actual 
programming  of  the  orbit  a  simple  matter  of  initializing  DISSPLA  for  the  type  of  mon- 
itor being  used,  setting  up  the  plotting  area,  initializing  the  axis  and  axis  scale,  and  then 
plotting  the  desired  curve  from  points  contained  in  arrays.  This  is  a  simplified  explana- 
tion of  DISSPLA,  but  for  further  details  on  DISSPLA  programming  refer  to  the 
DISSPLA  user's  manual  [Ref.  5].  DISSPLA  also  supplies  subroutines  to  draw  a  variety 
of  projections  of  the  world  and  fill  the  projections  with  coast  lines,  latitude  lines  and 
longitude  lines.  There  are  a  couple  of  DISSPLA  requirements  that  did  require  special 
handling  in  the  program.  The  requirement  that  the  data  be  supplied  in  arrays  forced  the 
program  to  load  arrays  with  the  required  position  and  parameters  and  to  keep  a  counter 
for  the  number  in  the  arrays.  The  array  format  requires  the  size  of  the  array  be  specified 
in  the  beginning  of  the  program.  The  array  size  needs  to  be  large  enough  to  hold  a 
number  of  orbits,  but  not  so  large  as  to  waist  storage  space.  The  program  will  continue 
to  add  orbital  data  to  the  arrays  until  the  user  chooses  to  delete  the  previous  orbits.  If 
a  new  initial  position  and  velocity  is  entered  or  if  the  arrays  will  overflow  with  the  next 
orbit  the  arrays  will  automatically  delete  all  previous  orbits.  DISSPLA  also  requires  that 
all  data  be  in  single  precision  format.  The  program  calculates  all  orbits  in  double  pre- 
cision in  order  to  limit  the  effect  of  round-off  error,  but  by  using  the  single  precision  data 
for  plotting  will  not  affect  the  accuracy  of  the  plot  in  any  way. 
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The  subroutines  used  to  display  the  orbits  will  be  covered  in  the  following  three 
sections; 

A.  PERIFOCAL  PLOT 

The  plotting  of  the  orbit  in  the  perifocal  coordinate  system  is  the  easier  of  the  two 
types  of  plots.  Since  the  perifocal  coordinate  system  has  the  orbital  plane  as  the  fun- 
damental plane,  the  only  requirements  to  describe  the  orbit  in  the  perifocal  coordinate 
system  are  arrays  with  the  true  anomaly  and  the  radius  to  the  satellite.  To  give  the  user 
a  sense  of  the  size  of  the  plot,  the  axis  length  varies  with  the  eccentricity  and  semi-major 
axis  length.  Also  a  plot  of  the  earth  is  plotted  to  the  same  scale,  with  the  pole  or  center 
of  the  plot  on  the  origin  of  the  axis.  The  latitude  of  the  earth  at  the  center  of  the  plot 
will  van."  with  the  inclination  of  the  orbit.  This  plot  will  allow  the  user  to  see  a  relative 
view  of  the  satellite's  coverage  in  the  minus  Z'  axis  direction  of  the  perifocal  coordinate 
system. 

B.  GROUND  TRACK 

The  ground  track  plot  is  a  very  complex  subroutine  compared  with  the  perifocal 
plot.  Because  the  ground  track  is  not  a  continuous  curve  a  procedure  to  handle  the 
satellite  ending  at  one  end  of  the  plot  and  wrapping  around  to  the  other  end  was  devel- 
oped. The  wrap  around  problem  is  avoided  in  most  orbits  by  plotting  the  orbit  in  seg- 
ments with  the  following  two  rules.  Each  segment  begins  at  the  beginning  of  a  new  plot 
or  at  the  edge  of  the  plot  area,  and  ending  when  the  satellite  would  wrap  around  to  the 
other  side  of  the  plot.  At  the  beginning  of  a  segment  if  the  position  of  the  satellite  is 
within  five  degrees  of  the  edge  of  the  plot,  that  position  and  any  other  positions  within 
that  five  degree  boundary  will  not  be  plotted.  The  segment  will  end  when  the  satellite 
is  within  ten  degrees  of  the  edge  of  the  plot.  The  above  restrictions  imposed  on  the 
segments  of  the  plot  will  not  substantially  affect  the  interpretation  or  usefulness  of  the 
plot.  The  ground  track  is  plotted  on  top  of  a  cylindrical  equidistant  projection  of  the 
world,  with  the  world  coast  lines  and  a  longitude-latitude  grid  for  reference. 

C.  DATA 

Information  concerning  the  orbit  is  displayed  on  the  lower  half  o[  the  plot.  The 
information  is  designed  to  supply  the  user  with  enough  of  the  basic  orbital  elements  and 
other  parameters  affecting  the  orbit  to  be  able  to  evaluate  what  basic  type  of  orbit  the 
satellite  is  in.  and  the  effects  of  velocity  changes  and  perturbing  forces  have  on  the  orbit. 
The  following  data  are  plotted:  inclination! i).  semi-major  axis  (a),  eccentricity  (e).  period 
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zee  and  perigee  velocity  and  radius,  average  time  rate-of-change  of  orbital  el- 
•    -■  average  magnitude  of  perturbing  forces  per  unit  mass. 


VII.     CONCLUSIONS  AND  RECOMMENDATIONS 

The  program  supplies  the  student  with  an  interactive  tool  to  study  the  orbital  mo- 
tion of  satellites  around  the  earth.  The  student  can  investigate  a  variety  of  orbits  by 
varying  the  orbital  parameters,  command  velocity  changes,  and  observe  the  effects  of 
perturbing  forces. 

The  student  is  provided  with  two  options  for  entering  the  initial  position  and  veloc- 
ity of  the  satellite.  The  program  could  be  expanded  to  provide  the  student  with  the  ad- 
ditional options  of  entering  either  orbital  parameters  or  a  ground  observation  data  and 
have  the  program  calculate  the  initial  position  and  velocity  from  this  data.  Also  the 
student  is  limited  to  orbits  with  eccentricities  less  than  one  (elliptic  orbits).  The  program 
could  be  also  be  expanded  to  include  more  eccentric  orbit  for  Lunar,  interplanetary,  and 
missile  trajectories.  The  perturbing  orbit  is  calculated  for  orbits  around  the  earth  with 
relatively  small  perturbing  forces  in  relation  to  the  earths  gravitational  force.  This  fact 
will  cause  the  program  to  produce  false  results  if  the  student  tries  to  calculate  lunar 
trajectories.  Special  routines  would  have  to  be  employed  when  the  perturbing  force  (the 
moons  gravitational  attraction)  is  comparable  to  the  earths  gravitational  attraction. 
This  will  not  become  a  factor  for  studying  current  satellite  orbits  out  to  the 
geosynchronous  radius  of  422-1 1.1  km. 

The  velocity  change  subroutines  move  the  satellite  to  a  location  close  to  the  desired 
location  before  a  velocity  change  is  imposed.  By  reducing  the  step  size  in  the  velocity 
change  subroutine,  this  error  could  be  reduced.  Precise  orbital  transfer  maneuvers  can 
be  modeled  by  reducing  this  error  caused  by  the  positioning  of  the  satellite  prior  to 
changing  the  velocity.  The  program  will  currently  provide  the  student  with  useful  plots 
for  gaining  experience  with  various  transfer  orbits  by  varying  the  magnitude  and  location 
of  the  velocity  changes. 

The  output  of  the  calculations  of  the  orbit  are  arrays  loaded  with  the  satellite's  po- 
sition and  select  orbital  parameters.  The  DISSPLA  subroutines  that  plot  the  points  are 
not  unique.  The  program  would  become  portable  to  personal  computers  with  these 
graphics  subroutines  written  in  FORTRAN  and  included  in  the  program. 

A  final  recommendation  is  that  the  display  of  the  ground  track  could  be  modified 
to  show  ground  coverage,  number  of  satellites  in  a  constellation,  and  other  elements 
necessary  for  planning  a  real-world  artificial  satellite  application. 
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APPENDIX  A.     ORBIT  PROGRAM 


PROGRAM  ORBIT 

THIS  PROGRAM  IS  AN  INTERACTIVE  TIME  STEP  SIMULATION  OF 
SATELLITES  AROUND  THE  EARTH.   PERTURBED  AND  UNPERTURBED  ORBITS 
ARE  CALCULATED  AND  PLOTTED.   VELOCITY  CHANGES  ARE  ALSO  PERMITTED 
AT  SPECIFIED  TRUE  ANOMALIES. 
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A  LIST 

A 

AL 

AP 

CHTA  = 

DT 

E 

EA 

EI 

EJ 

EK 

FR 

FS    = 

FW 

H 

HI 

HJ 

HK 

I 

IOPTl= 

IOPT2= 

CHANGE 

LAN  = 

LP 

MA 

MM 

MU 

N 

NI 

NJ 

NK 

NUM  = 

P 

PER  = 

PI 

RA 

RE 

R 

RI 

RJ 

RK 

T 

TA 

TDA  = 

TDAP  = 


OF  VARIABLES  USED  BY  THE  MAIN  PROGRAM  FOLLOWS: 

SEMI -MAJOR  AXIS 

ARGUMENT  OF  LONGITUDE 

ARGUMENT  OF  PERIGEE 

VELOCITY  CHANGE  LOCATION  TRUE  ANOMALY 

TIME  STEP 

ECCENTRICITY 

ECCENTRIC  ANOMALY 

I  VECTOR  OF  ECCENTRICITY 

J  VECTOR  OF  ECCENTRICITY 

K  VECTOR  OF  ECCENTRICITY 

R  VECTOR  OF  TOTAL  FORCE 

S  VECTOR  OF  TOTAL  FORCE 

W  VECTOR  OF  TOTAL  FORCE 

ANGULAR  MOMENTUM 

I  VECTOR  OF  ANGULAR  MOMENTUM 

J  VECTOR  OF  ANGULAR  MOMENTUM 

K  VECTOR  OF  ANGULAR  MOMENTUM 

INCLINATION 

PERTURBED  OR  UNPERTURBED  OPTION 

OPTIONS:  PLOT  NEXT  ORBIT,  CHANGE  INITIAL  VALUES, 

VELOCITY,  PLOT  ANOTHER  VIEW  OF  ORBIT,  QUIT 

LONGITUDE  OF  ASCENDING  NODE 

LONGITUDE  OF  PERIGEE 

MEAN  ANOMALY 

MEAN  MOTION 

GRAVITATIONAL  PARAMETER 

ASCENDING  NODE 

I  VECTOR  OF  ASCENDING  NODE 

J  VECTOR  OF  ASCENDING  NODE 

K  VECTOR  OF  ASCENDING  NODE 

STEP  COUNTER 

SEMI-LATUS  RECTUM 

PERIOD  OF  ORBIT 

PI 

RADIUS  OF  APOGEE 

RADIUS  OF  EARTH 

ORBITAL  RADIUS 

I  VECTOR  OF  ORBITAL  RADIUS 

J  VECTOR  OF  ORBITAL  RADIUS 

K  VECTOR  OF  ORBITAL  RADIUS 

TIME  COUNTER  IN  ORBIT 

TRUE  ANOMALY 

TOTAL  CHANGE  IN  SEMI -MAJOR  AXIS 

TOTAL  CHANGE  IN  ARGUMENT  OF  PERIGEE 


ORBOOOIO 
ORB00020 
ORB00030 
ORB00040 
ORB00050 
ORB00060 
ORB00O70 
ORB00080 
ORB00090 
ORBOOIOO 
ORB00110 
ORB00120 
ORB00130 
ORB00140 
ORB00150 
ORB00160 
ORB0017O 
ORB00180 
ORB00190 
ORB00200 
ORB00210 
ORB00220 
ORB00230 
0RB00240 
ORB00250 
ORB00260 
ORB00270 
0RB00280 
ORB00290 
ORB00300 
ORB00310 
ORB00320 
ORB00330 
ORB00340 
ORB00350 
ORB00360 
ORB00370 
ORBOO380 
ORB00390 
ORB00400 
ORB00410 
ORB00420 
0RB00430 
ORB00440 
0RB00450 
ORB00460 
ORB00470 
ORB00480 
ORB00490 
ORB005O0 
ORB00510 
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TDE  =  TOTAL  CHANGE  IN  ECCENTRICITY 

TDK  =  TOTAL  CHANGE  IN  ANGULAR  MOMENTUM 

TDI   =  TOTAL  CHANGE  IN  INCLINATION 

TDMA  =  TOTAL  CHANGE  IN  MEAN  ANOMALY 

TDMM  =  TOTAL  CHANGE  IN  MEAN  MOTION 

TDLAN=  TOTAL  CHANGE  IN  LONGITUDE  OF  ASCENDING  NODE 

TF   =  TIME  OF  FLIGHT 

TFDRA=  TOTAL  FORCE  OF  DRAG 

TFEA  =  TOTAL  FORCE  OF  EARTH'S  OBLATENESS 

TFMO  =  TOTAL  FORCE  FROM  MOON 

TFSU  =  TOTAL  FORCE  FROM  SUN 

TL   =  TRUE  Longitude  AT  EPOCH 

TT   =  TRUE  TIME  SINCE  SATELLITE  HAS  BEEN  IN  ORBIT 

V  =  SATELLITE  VELOCITY 

VI  =  I  VECTOR  OF  SATELLITE  VELOCITY 
VJ  =  J  VECTOR  OF  SATELLITE  VELOCITY 
VK   =  K  VECTOR  OF  SATELLITE  VELOCITY 

A  LIST  OF  THE  ARRAYS  USED  FOLLOWS: 

AINRAY  =  INCLINATION 

APRAY  =  ARGUMENT  OF  PERIGEE 

RARAY  =  RADIUS 

RIRAY  =  I  VECTOR  OF  RADIUS 

RJRAY  =  J  VECTOR  OF  RADIUS 

RKRAY  =  K  VECTOR  OF  RADIUS 

TARAY  =  TRUE  ANOMALY 

TIMRAY  =  TIME 

A  LIST  OF  SUBROUTINES  CALLED  BY  THE  MAIN  PROGRAM  WILL  FOLLOW: 

CALCEL  =  CALCULATES  THE  ORBITAL  ELEMENTS 

CHGVEL  =  ALLOW  THE  USER  TO  CHANGE  THE  VELOCITY  OF  THE  SATELLITE 

INPUTS  =  PROMPTS  USER  FOR  INITIAL  POSITION  AND  VELOCITY 

INTSUM  =  INITIALIZES  THE  SUMS  IN  THE  ARRAYS 

NEWELT  =  CALCULATE  NEW  ORBITAL  ELEMENTS  FROM  TIME  STEP 

NEWPOS  =  CALCULATE  NEW  POSITION  VECTOR 

NEWVEL  =  CALCULATE  NEW  VELOCITY  VECTOR 

OPTION  =  GIVE  THE  USER  THE  OPTIONS  Permitted  IN  THE  PROGRAM 

PLOTS  =  PLOTS  THE  ORBITS 

PRETUR  =  CALCULATES  THE  PERTURBED  ORBIT 

STORE  =  STORE  THE  POSITION  DATA  IN  ARRAYS 

UNPRET  =  CALCULATE  THE  UNPERTURBED  ORBIT 

BEGIN  MAIN  PROGRAM 

DOUBLE  PRECISION  PI ,MU,RI ,RJ,RK,R,VI , VJ,VK,V,HI ,HJ,HK,H, 
+  NI.NJ.NK.N.P.EI.EJ.EK.E.A.I.LAN.AP.TA.AL.LP.TL.PER.EA, 
+   MM,MA,T,DT,TF,  FR,FS,FW,TT,CHTA,RA,VA,TEMPTA,RE 

DIMENSION  TARAY(500),RARAY(500),RIRAY(500),RJRAY(500),RKRAY(500)5 
+   AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* 1 , LOOP , YORN , ORLOOP 

PI  =  3.  141592653589794 


ORB00520 
ORB00530 
0RB00540 
ORB00550 
ORB00560 
ORB00570 
ORBO0580 
ORB0059O 
ORB0O600 
ORB00610 
ORB00620 
ORB00630 
ORB00640 
ORB00650 
ORB00660 
ORBO067O 
ORB00680 
ORB00690 
ORB00700 
ORB00710 
ORB0072O 
ORB00730 
ORB00740 
ORB00750 
ORB00760 
ORB00770 
ORB00780 
ORB00790 
ORB00800 
ORB00810 
ORB00820 
ORB00830 
ORB00840 
ORB00850 
ORB00860 
ORB00870 
ORB00880 
ORB00890 
ORB00900 
ORB00910 
ORB00920 
ORB00930 
ORB00940 
ORB0O950 
ORB00960 
ORB0097O 
ORB00980 
ORB00990 
ORBOIOOO 
ORB0101O 
ORB01020 
ORB01030 
ORB01040 
ORB01050 
ORB01060 
ORB01070 
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MU  =  3. 986012D+05 
RE  =  6. 378145D+03 

USER  INTRO  TO  PROGRAM 
CALL  INTRO 

■     ENTERED  MAIN  PROGRAM  LOOP 

LOOP  =  'Y' 
10    IF  (LOOP  .  EQ.  'Y')  THEN 

INITIALIZE  STEP  COUNTER  AND  TRUE  TIME 
20      NUM  =  1 
TT  =  0.  0 

PROMPT  USER  FOR  INITIAL  POSITION  AND  VELOCITY 
CALL  INPUTS(RI,RJ,RK,R,VI,VJ,VK,V,MU,LOOP,PI) 

EXIT  PROGRAM 

IF  (LOOP  .EQ.  'N')  THEN 

GOTO  10 
END  IF 

•        CALCULATE  AND  STORE  ORBITAL  ELEMENTS 

CALL  CALCEL(RI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E,A,I,LAN, 
+  LP,TA,PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

CALL  STORE (RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY,TARAY, 
+  NUM,I,AP,AINRAY,APRAY,TT,TIMRAY) 

PRINT  DATE  FOR  USER  TO  REVIEW 


PRINT- 
PRINT* 
PRINTS- 
PRINT* 

PRINTS- 
PRINT^ 
PRINT* 
PRINT* 
PRINT* 
DEGI  = 
PRINT* 
PERHRS 
PRINT* 
PRINT* 


VI 
VJ 
VK 

V 
RI 
RJ 
RK 

R 


VI, 
VJ. 
VK, 
V,' 
RI, 
RJ, 
RK, 
R,' 


KM/S' 

KM/S' 

KM/S' 
KM/S' 

KM' 

KM' 

KM* 
KM' 


DEGREES' 


ENTER 


"ylt 


OR 


"N" 


50 


'ECCENTRICITY  =' ,E 
SNGL((180.0/PI)*I) 
'INCLINATION  =' ,DEGI, 
=  SNGL( PER/3600.  0) 
'PERIOD  =' ,PERHRS,'  HOURS' 
'ARE  THESE  VALUES  CORRECT? 

READ''-  f  yorn 

CALL  EXCMS('CLRSCRN') 

IF  (.NOT.  YORN  .  EQ.  'Y')  THEN 
GOTO  20 

END  IF 


CALCULATE  TIME  STEP  AND  SET  TIMER  TO  ONE  TIME  STEP 
DT  =  PER/50 
T  =  DT 

STEP  SATELLITE  TO  PERIGEE  POINT  AND  RECORD 
IF  ((TA.  GT.  0.  063).  AND.  (TA.  LT.  6.  21))  THEN 


TT  = 


+  DT 


ORB01080 
ORB01090 
ORB01100 
ORB01110 
ORB01120 
ORB01130 
ORB01140 
ORB01150 
ORB01160 
ORBO1170 
ORB01180 
ORB01190 
ORB01200 
0RB01210 
ORB01220 
ORB01230 
ORB01240 
ORB01250 
ORB01260 
ORB01270 
ORB01280 
ORB01290 
ORB01300 
ORB01310 
ORB01320 
ORB01330 
ORB01340 
ORB01350 
ORB01360 
ORB01370 
ORB01380 
ORB01390 
ORB01400 
0RB01410 
0RB01420 
0RB01430 
ORB01440 
ORB01450 
0RB01460 
ORB01470 
ORB01480 
0RB01490 
ORB01500 
0RB01510 
ORB01520 
0RB01530 
ORB01540 
ORB01550 
ORB01560 
ORB01570 
ORB01580 
ORB01590 
ORB01600 
ORB01610 
ORB01620 
ORB01630 
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CALL  NEWELT(MM,MA,E,EA,TA,TF,DT,PI,PER) 
CALL  NPOS(RI,RJ,RK,R,LAN,AP,I,  TA,A,E) 
CALL  NVEL( E , P , TA , LAN , AP , I , VI , VJ , VK , V , MU ) 
NUM  =  NUM  +  1 

CALL  STORE( RI , RJ , RK , R , TA , RIRAY , R JRAY , RKRAY , RARAY , TARAY , 
+  NUM ,  I ,  AP , AINRAY , APRAY , TT , TIMRAY) 

T  =  T  +  DT 
GOTO  50 
END  IF 

'  CALCULATE  ELEMENTS  FROM  PERIGEE  POINT 

CALL  CALCEL(RI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E,A,I,LAN, 
+  LP,TA,PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

DT  =  PER/50 

T  =  DT 

STORE  FIRST  Unperturbed  ORBIT 

CALL  UNPRET( DT , PER , AL , LAN , AP , I , RI , R J , RK , R , VI , VJ , VK , V , 
+  MU,PI,H,A,E,N,TA,P,MM,MA,EA,TF,T,NUM,RIRAY,RJRAY, 

+  RKRAY, RARAY, TARAY, AINRAY, APRAY , TIMRAY, TT) 

■  INITIALIZE  SUMS  FOR  FORCE  AND  ORBITAL  ELEMENT  CHANGES  TO  ZERO 
CALL  INTSUM( TFEA , TFSU , TFMO , TFDRA , TDI , TDA , TDE , TDMM , TDMA , TDLAN , 

+  TDH,TDAP) 

PLOT  FIRST  UNPERTURBED  ORBIT 

7 0      CALL  PLOTS (RIRAY , R JRAY , RKRAY , RARAY , TARAY , NUM , P I , I , LP , A , E , TF , 
+  AINRAY, APRAY, TIMRAY, TFEA, TFSU, TFMO, TFDRA, PER, TDI, TDA 

+  TDE , TDMM , TDMA , TDLAN , TDH , TDAP , MM , MA , LAN , H , AP , R , V ) 

BEGIN  NEW  ORBIT  OPTIONS 
IOFT1  =  1.  Unperturbed  ORBIT 

=  2.  Perturbed  ORBIT 

=  3.  QUIT 

■  I0PT2  =  1.  PLOT  NEXT  ORBIT 

■  =  2.  CHANGE  INITIAL  VALUES 

•  =  3.  CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  Anomaly 

•  =  4.  PLOT  ANOTHER  VIEW  OF  SAME  ORBIT 

ALSO  ASKED  IF  WANT  TO  CLEAR  ALL  PREVIOUS  ORBITS 

CALCULATE  ELEMENTS  AT  PERIGEE 
80      CALL  CALCEL(RI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E,A,I,LAN, 

+  LP,TA,PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

•  CHECK  FOR  POSSIBLE  ARRAY  OVERFLOW 
IF  (NUM  . GT.  425)  THEN 

PRINT-"-, 'ARRAYS  ARE  FULL' 

PRINT*,' PREVIOUS  ORBITS  WILL  BE  ERASED!  ' 
NUM  =  1 

CALL  STORE ( R I, R J, RK, R, TA, R I RAY,R JRAY, RKRAY, RARAY, TARAY, 
+  NUM, I, AP, AINRAY, APRAY, TT, TIMRAY) 

END  IF 

PROMPT  USER  FOR  DESIRED  OPTION 


ORB01640 
ORB01650 
ORB01660 
ORB01670 
0RB01680 
ORB01690 
ORB01700 
ORB01710 
ORB01720 
ORB01730 
ORB01740 
ORB01750 
ORB01760 
ORB01770 
ORB01780 
ORB01790 
ORB01800 
ORB01810 
ORB01820 
ORB01830 
ORB01840 
ORB01850 
ORB01860 
ORB01870 
ORB01880 
ORB01890 
ORB01900 
ORB01910 
,ORB01920 
ORB01930 
0RB01940 
ORB01950 
ORB01960 
ORB01970 
ORB01980 
ORB01990 
ORB02000 
ORB02010 
ORB02020 
ORB02030 
ORB02040 
ORB02050 
ORB02060 
ORB02070 
ORB02080 
ORB02090 
ORB02100 
ORB02110 
ORB02120 
ORB02130 
0RB02140 
ORB02150 
ORB02160 
ORB02170 
ORB02180 
ORB02190 


* 


* 
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CALL  0PTI0N(I0PT1,I0PT2, NUM, RIRAY, RJRAY, RKRAY, RARAY, 
+  TARAY , AINRAY , APRAY , TIMRAY ) 

Initialize  SUMS  FOR  FORCE  AND  ORBITAL  ELEMENT  CHANGES  TO  ZERO 
CALL  INTSUM( TFEA , TFSU , TFMO , TFDRA , TDI , TDA , TDE , TDMM ,TDMA , TDLAN , 
+  TDH,TDAP) 

SET  TIME  COUNTER  TO  ONE  TIME  STEP 
T  =  DT 

OPTION:  PLOT  THE  NEXT  ORBIT 
IF  (I0PT2  .EQ.  1)  THEN 

CALCULATE  AND  PLOT  UNPERTURBED  ORBIT 

IF(IOPTl  .EQ.  1)  THEN 

CALL  UNPRETC DT , PER , AL , LAN , AP , I , RI , R J , 
+  RK,R,VI,VJ,VK,V,MU,PI,H,A, 

+  E , N , TA , P , MM , MA , E A , TF , T , NUM , R I RAY , R JRAY , RKRAY , 

+  RARAY, TARAY, AINRAY, APRAY, TIMRAY, TT) 

CALL  PLOTS ( RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 
+  PI, I, LP, A, E,TF, AINRAY, APRAY, TIMRAY, 

+  TFEA, TFSU, TFMO, TFDRA, PER, 

+  TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

+  MM,MA,LAN,H,AP,R,V) 

CALCULATE  AND  PLOT  PERTURBED  ORBIT 

ELSEIF(IOPTl  .EQ.  2)  THEN 

CALL  PRETURC  DT , PER , AL , LAN , AP , I , 
+  RI,RJ,RK,R,VI,VJ,VK,V,FR,FS,FW, 

+  M U ,  P I ,  H ,  A ,  E  ,  N ,  T A ,  P ,  MM ,  M A ,  E  A ,  TF ,  T ,  N UM , 

+  R IRAY , RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , 

+  TIMRAY, TT, TFEA, TFSU, TFMO, TFDRA, 

+  TDI , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP ) 

CALL  PLOTS (RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 
+  P I, I, LP, A, E,TF, AINRAY, APRAY, TIMRAY, 

+  TFEA, TFSU, TFMO, TFDRA, PER, 

+  TD I , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

+  MM,MA,LAN,H,AP,R,V) 

END  IF 

GOTO  THE  BEGINNING  OF  THE  PROGRAM  TO  CHANGE  THE  INITIAL  VALUES 
ELSEIF  (IOPT2  . EQ.  2)  THEN 
GOTO  20 

CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  ANOMALY  AND 

PLOT  THE  NEW  ORBIT 

ELSEIF  (IOPT2  . EQ.  3)  THEN 

CALL  CHGVEL(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R, 
+  VI,VJ,VK,V,MU,PI, 

+  H,A,E,N,TA,P,MM,MA,EA,TF,T,NUM,RIRAY, 

+  RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , 

+  TIMRAY,TT,EI,EJ,EK,LP,HI,HJ,I0PT1, 

+  TFEA, TFSU, TFMO, TFDRA, TDI, TDA, TDE, TDMM, 

+  TDMA, TDLAN, TDH, TDAP) 

CALL  PLOTS ( RIRAY , RJRAY , RKRAY , RARAY , TARAY , NUM , 
+  PI , I , LP , A, E ,TF , AINRAY , APRAY , TIMRAY , 


0RB02200 
0RB02210 
0RB02220 
ORB02230 
0RB02240 
ORB02250 
ORB02260 
ORB02270 
ORB02280 
ORB02290 
ORB02300 
ORB02310 
ORB02320 
ORB02330 
0RB02340 
ORB02350 
ORB02360 
ORB02370 
ORB02380 
ORB02390 
ORB02400 
0RB02410 
ORB02420 
ORB02430 
0RB02440 
ORB02450 
ORB02460 
ORB02470 
0RB02480 
0RB02490 
ORB02500 
ORB0251O 
0RB02520 
ORB02530 
ORB02540 
ORB02550 
ORB 025 60 
ORB02570 
ORB02580 
ORB02590 
ORB02600 
ORB02610 
ORB02620 
ORB02630 
ORB02640 
ORB02650 
ORB02660 
ORB02670 
ORB02680 
ORB02690 
ORB02700 
ORB02710 
ORB02720 
ORB02730 
ORB02740 
ORB02750 
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+ 
+ 
+ 


+ 
+ 
+ 

+ 


90 


TFE A , TFSU , TFMO , TFDRA , PER , 

TD I , TD A , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

MM,MA,LAN,H,AP,R,V) 

PLOT  ANOTHER  VIEW  OF  THE  SAME  ORBIT 
ELSEIF  (IOPT2  .  EQ.  4)  THEN 

CALL  PLOTS (RIRAY,RJRAY,RKRAY,RARAY,TARAY,NUM, 

PI,I,LP,A,E,TF,AINRAY,APRAY,TIMRAY, 

TFEA , TFSU , TFMO , TFDRA , PER , 

TD I , TDA , TDE , TDMM , TDMA , TDLAN , TDH , TDAP , 

MM,MA,LAN,H,AP) 

STOP  THE  PROGRAM 

ELSEIF  (IOPT2  . EQ.  5)  THEN 

GOTO  90 
ELSE 

'  PR I NT*,* INVALID  ENTRY!  ' 

GOTO  80 
END  IF 

CHECK  IF  SATELLITE  Impacted  THE  EARTH  AND  GO  TO  THE  BEGINNING 

IF  (R  .LE.  6450.0)  THEN 

PRINT*, 'SATELLITE  WILL  IMPACT  THE  EARTH! !  !  ' 
PR I NT*,' PROGRAM  WILL  RESET  TO  THE  BEGINNING' 
GOTO  20 

END  IF 

GOTO  THE  TOP  OF  THE  OPTION  LOOP 
GOTO  80 

GIVE  THE  USER  A  CHANCE  TO  RECOVER  THE  PROGRAM 

PRINT"-, 'THIS  IS  YOUR  LAST  CHANCE!  ' 

PRINT-, 'DO  YOU  WANT  TO  CONTINUE?' 

PRINT--", 'AND  GOTO  THE  Beginning  OF  THE  PROGRAM?' 


PRINT--', 'ENTER  "Y   OR  "N"  : 
READ'--,  LOOP 
PRINT* , LOOP 
GOTO  10 
END  IF 

*  DISSPLA  SUBROUTINE  TO  TELL  GRAPHICS  TERMINAL  PLOTTING 
SESSION  IS  DONE 

CALL  DONEPL 

STOP 

END 

SUBROUTINE  INTRO 

*  THIS  SUBROUTINE  WILL  GIVE  THE  USER  A  Brief  INTRODUCTION  OF  THE 

USES  OF  THE  PROGRAM 

PRINT-', 'THIS  PROGRAM  IS  A  GRAPHICS  DISPLAY  OF  Satellite  ORBITS. 
PRINT--", 'YOU  WILL  BE  ASKED  TO  INPUT  THE  INITIAL  VELOCITY  AND' 
PRINT---, 'POSITION  VECTORS  OF  THE  Satellite.   THE  PROGRAM  WILL   ' 
PRINT'--, 'THEN  CALCULATE  THE  ORBITAL  PARAMETERS  AND  THE  ' 


ORB02760 
ORB02770 
0RB02780 
ORB02790 
ORB02800 
ORB02810 
ORBO2820 
ORB02830 
0RB02840 
ORB02850 
ORB02860 
ORB02870 
ORB02880 
ORB02890 
ORB02900 
ORB02910 
ORB02920 
ORB02930 
0RB02940 
ORB02950 
ORB02960 
ORB02970 
ORB02980 
ORB02990 
ORB03000 
ORB03010 
ORB03020 
ORB03030 
ORB03040 
ORB03050 
ORB03060 
ORB03070 
ORB03080 
ORB03090 
ORB03100 
ORB03110 
ORB03120 
ORB03130 
0RB03140 
ORB03150 
ORB03160 
ORB03170 
ORB03180 
ORB03190 
ORB03200 
ORB03210 
ORB03220 
ORB03230 
ORB03240 
ORB03250 
ORB03260 
ORB03270 
ORB03280 
ORB03290 
ORB03300 
ORB03310 
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PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT'"' 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
RETURN 
END 


Unperturbed  ORBIT.   THE  USER  WILL  THEN  HAVE  THE' 

CHOICE  OF  DISPLAYS: ' 

-PERIFOCAL  (SHOWS  RELATIVE  SIZE  OF  ORBIT)' 
-Equatorial  (SHOWS  ORBIT  INCLINED,  USER  INPUT' 
LONGITUDE  TO  VIEW  AT)' 

-GROUND  TRACK' 

t 

THE  USER  IS  THEN  ASKED  TO  CHOOSE  ONE  OF  THE  FOLLOWING: 

-Unperturbed  ORBITS' 

-Perturbed  ORBITS' 

-VELOCITY  CHANGES' 
THE  USER"S  CHOICE  WILL  BE  USED  IN  DEVELOPING  THE' 

GRAPHICAL  OUTPUT. ' 

i 

THE  USER  IS  THEN  GIVEN  THE  FOLLOWING  CHOICES: ' 
-CLEAR  ALL  THE  PREVIOUS  ORBITS' 
-CHANGE  THE  INITIAL  PARAMETERS' 
-CHANGE  VELOCITY  AT  A  SPECIFIC  TRUE  Anomaly' 
-PLOT  ANOTHER  VIEW  OF  THE  SAME  ORBIT' 


:Wr*****Vry«V}Vyr***********iWf  ************** 

SUBROUTINE  OPTION( IOPT1 , IOPT2 ,NUM,RIRAY,RJRAY,RKRAY,RARAY, 
+  TARAY,AINRAY,APRAY,TIMRAY) 

THIS  SUBROUTINE  GIVES  THE  USER  A  CHOICE  OF  OPERATIONS  THAT  CAN  BE 
PERFORMED  ON  THE  PROGRAM  AND  RETURNS  THE  USERS  CHOICE  WITH 
VARIABLES  IOPT1  AND  IOPT2 

DIMENSION  RIRAY(500),RJRAY(500),RKRAY(500),RARAY(500),TARAY(500), 
+   AINRAY(500),APRAY(500),TIMRAY(500) 
CHARACTER* l,YORN 
IOPT1  =  0 

PROMPT  USER  FOR  OPTION 
103   PRINT*, 'WHICH  OF  THE  FOLLOWING  OPTIONS  WOULD  YOU  LIKE: ' 
PRINT*,'   1. -CALCULATE  THE  NEXT  ORBIT  USING  THE  SAME* 
PRINT*,'      PARAMETERS' 

PRINT*,'   2. -CHANGE  THE  INITIAL  PARAMETERS  OF  THE  ORBIT' 
PRINT*,'   3. -CHANGE  THE  VELOCITY  AT  A  POINT  IN  THE  ORBIT' 
PRINT*,'      (THE  UNPERTURBED  ORBIT  WILL  BE  USED)' 
PRINT*,'   4. -PLOT  ANOTHER  VIEW  OF  THE  ORBIT(S)' 
PRINT'S'   5. -QUIT' 
PRINT*, 'ENTER  1,  2,  3,  4,  OR  5: ' 
READ*,I0PT2 
PRINT*, IOPT2 
CALL  EXCMS( 'CLRSCRN') 
IF  (  IOPT2  .GT.  5)  THEN 

GOTO  103 
END  IF 

*     Prompt  USER  FOR  TYPE  OF  ORBIT  DESIRED 
105   IF  (I0PT2  .  EQ.  1)  THEN 

PRINT*,' WHICH  TYPE  OF  ORBIT  WOULD  YOU  LIKE  TO  SEE,' 
PRINT*,'    1. -Unperturbed  ORBITS' 


ORB03320 
ORB03330 
0RB03340 
ORB03350 
ORB03360 
ORB03370 
ORB03380 
ORB03390 
ORB03400 
ORB03410 
ORB03420 
ORB03430 
ORB03440 
ORB03450 
0RB03460 
ORB03470 
ORB03480 
ORB03490 
ORB03500 
ORB03510 
ORB03520 
ORB03530 
ORB03540 
ORB03550 
ORB03560 
ORB03570 
ORB03580 
ORB03590 
ORB03600 
ORB03610 
ORB03620 
ORB03630 
ORB03640 
ORB03650 
ORB03660 
ORB03670 
ORB03680 
ORB03690 
ORB03700 
ORB03710 
ORB03720 
ORB03730 
ORB03740 
ORB03750 
ORB03760 
ORB03770 
ORB03780 
ORB03790 
ORB03800 
ORB03810 
ORB03820 
ORB03830 
ORB03840 
ORB03850 
0RB03860 
ORB03870 


PRINT*,1    2.  -Perturbed  ORBITS' 

PRINT*,'  ENTER  1  OR  2: ' 

READ*, I OPT 1 

PRINT*, IOPT1 

CALL  EXCMS('CLRSCRN') 

IF  ((IOPT1  . NE.  1)  .AND.  (IOPT1  .  NE.  2))  THEN 

PRINT*, 'INVALID  ENTRY! ' 

GOTO  105 
END  IF 
END  IF 

PROMPT  USER  TO  CLEAR  PREVIOUS  ORBITS 
107   IF  ((IOPT2  . EQ.  1)  .OR.  (IOPT2  .  EQ.  3))  THEN 

PRINT*,' DO  YOU  WANT  TO  CLEAR  THE  PREVIOUS  ORBITS? 

PRINT*, 'ENTER   "Y"   OR  "N"    : ' 

READ*,YORN 

PRINT'S  YORN 

CALL  EXCMS( 'Clrscrn') 

IF   (YORN   .EQ.     'Y')   THEN 

RIRAY(l)  =  RIRAY(NUM) 
RJRAY(NUM) 
RKRAY(NUM) 
RARAY(NUM) 
TARAY(NUM) 
=  AINRAY(NUM) 


RJRAY(l)  : 
RKRAY(l)  = 
RARAY(l)  : 
TARAY(l)  : 
AINRAY(l) 


BE  CAPITOL  LETTERS 


NE. 
.NE 


2).  AND.  (I0PT2 
5))  THEN 


NE.  3)  .AND. 


APRAY(l)  =  APRAY(NUM) 
TIMRAY(l)  =  TIMRAY(NUM) 
NUM  =  1 

ELSEIF  (YORN  .  NE.  'N')  THEN 
PRINT*, ' INVALID  ENTRY!  !  ' 
PRINT*, 'ALL  INPUTS  MUST 
GOTO  107 

END  IF 
END  IF 

CHECK  FOR  INVALID  OPTION 

IF  ((IOPT2  .NE.  1).AND.  ( IOPT2  . 
+     (IOPT2  .NE.  4).  AND.  ( IOPT2 
PRINT*, 'INVALID  ENTRY! ' 
GOTO  103 
END  IF 
RETURN 
END 


*  COORDINATE  TRANSFORMATIONS 

yoYyr-yiYV—Vy-y-iYycVryryrVrycVfyryTyryryriY^^ 

SUBROUTINE  PQWIJK(LAN, AP, INC,P,Q,W, I , J,K) 

*  THIS  SUBROUTINE  TRANSFORMS  PQW  COORDINATES  TO  UK  COORDINATES 

DOUBLE  PRECISION   INC,P,Q,W, I , J,K,R11 ,R12 ,R13 ,R21 ,R22 ,R23, 
+    R31,R32,R33,LAN,AP 
Rll  =  DCOS(LAN)*DCOS(AP)  -  DSIN(LAN)*DSIN( AP)*DCOS(  INC") 
R12  =  -DCOS(LAN)*DSIN(AP)  -  DSIN( LAN)*DCOS( AP)*DCOS( INC) 
R13  =  DSIN(LAN)*DSIN(INC) 


ORB03880 
0RB03890 
ORB03900 
ORB03910 
0RB03920 
ORB03930 
ORB03940 
ORB03950 
ORB03960 
ORB03970 
ORB03980 
ORB03990 
ORB04000 
ORB04010 
ORB04020 
ORB04030 
ORB04040 
ORB04050 
ORB04060 
ORB04070 
ORB04080 
ORB04090 
ORB04100 
0RB04110 
0RB04120 
ORB04130 
ORB04140 
0RB04150 
0RB04160 
ORB04170 
0RB04180 
ORB04190 
0RB04200 
ORB04210 
0RB04220 
ORB04230 
ORB04240 
0RB04250 
ORB04260 
ORB04270 
0RB04280 
ORB04290 
ORB04300 
ORB04310 
ORB04320 
ORB04330 
ORB04340 
ORB04350 
ORB04360 
ORB04370 
ORB04380 
ORB04390 
ORB04400 
ORB04410 
ORB04420 
ORB04430 
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DSIN(LAN)*DCOS(AP)  +  DCOS(LAN)*DSIN(AP)*DCOS( INC) 

-DSIN(LAN)*DSIN(AP)  +  DCOS(LAN)*DCOS(AP)*DCOS( INC) 

-DCOS(LAN)*DSIN(INC) 

DSIN(AP)*DSIN(INC) 

DCOS(AP)*DSIN(INC) 

DCOS(INC) 

1*P  +  R12*Q  +  R13*W 

1*P  +  R22*Q  +  R23*W 

1*P  +  R32*Q  +  R33*W 


SUBROUTINE  IJKPQW(LAN,AP, INC , I , J,K,P,Q ,W) 
*     THIS  SUBROUTINE  TRANSFORMS  IJK  COORDINATES  TO  PQW  COORDINATES 


R21 

=  : 

R22 

= 

R23 

= 

R31 

=  ; 

R32 

R33 

I  = 

Rl 

J  = 

R2 

K  = 

R3 

RETURN 

END 

DOU 
h 

Rll 
R21 
R31 
R12 
R22 
R32 
R13 
R23 
R33 
P  = 

Q  = 
W  = 
RET 
END 


BLE  PRECISION   INC , I , J,K,P,Q,W,R11 ,R12 ,R13 ,R21 ,R22 ,R23 , 
R31,R32,R33,LAN,AP 

=  DCOS(LAN)*DCOS(AP)  -  DSIN(LAN)*DSIN( AP)*DCOS( INC) 

=  -DCOS(LAN)*DSIN(AP)  -  DSIN(LAN)*DCOS(  AP)*DCOS(  INC) 

=  DSIN(LAN)*DSIN(INC) 

=  DSIN(LAN)*DCOS(AP)  +  DCOS( LAN)*DSIN( AP)*DCOS( INC) 

=  -DSIN(LAN)--'-'DSIN(AP)  +  DCOS(LAN)*DCOS( AP)*DCOS(  INC) 

=  -DCOS(LAN)*DSIN(INC) 

=  BSIN(AP)*DSIN(INC) 

=  DCOSCAP)*DSIN(INC) 

=  DCOS(INC) 

R11*I  +  R12'-'-J  +  R13*K 

R21-I  +  R22-J  +  R23-''K 

R31*I  +  R32*J  +  R33*K 
URN 


SUBROUTINE  IJKRSW( LAN, AL, INC , I , J,K,R,S ,W) 

THIS  SUBROUTINE  CHANGES  FROM  IJK  COORDINATES  TO  RSW  COORDINATES 


DOUBLE 
+   R31 
Rll  = 
R12  = 
R13  = 
R21  = 
R22  = 
R23  = 
R31  = 
R32  = 
R33  = 
R  =  Rl 
S  =  R2 
W  =  R3 
RETURN 
END 


PRECISION   INC,I,J,K,R,S,W,R11,R12,R13,R21,R22,R23, 
,R32,R33,LAN,AL 

DCOS(LAN)*DCOS(AL)  -  DSIN( LAN)*DCOS( INC)*DSIN( AL) 
DSIN(LAN)*DCOS(AL)  +  DSIN(AL)*DCOS(LAN)*DCOS(INC) 
DSIN(INC)*DSIN(AL) 

-DCOS(LAN)*DSIN(AL)-DSIN(LAN)*DCOS(INC)*DCOS(AL) 
-DSIN(LAN)*-'DSIN(AL)  +  DCOS(LAiN)*DCOS(INC)*DCOS(AL) 
DSIN(INC)*DCOS(AL) 
DSIN(LAN)*DSIN(INC) 
-DCOS(LAN)*DSIN(INC) 
DCOS(INC) 

1*1  +  R12*J  +  R13*K 
1*1  +  R22*J  +  R23*K 
1*1   +  R32*J  +  R33*K 


0RB04440 
0RB04450 
ORB04460 
ORB04470 
0RB04480 
0RB04490 
0RB04500 
0RB04510 
ORB04520 
ORB04530 
ORB04540 
ORB04550 
ORB04560 
ORB04570 
0RB04580 
ORB04590 
0RB04600 
0RB04610 
0RB04620 
ORB04630 
ORB04640 
0RB04650 
ORB04660 
ORB04670 
ORB04680 
ORB04690 
ORB04700 
ORB04710 
ORB04720 
ORB04730 
ORB04740 
ORB04750 
ORB04760 
ORB04770 
ORB04780 
ORB04790 
ORB04800 
ORB04810 
ORB04820 
ORB04830 
0RB04840 
ORB04850 
ORB04860 
ORB04870 
ORB04880 
ORB04890 
ORB04900 
ORB04910 
ORB04920 
ORB04930 
ORB04940 
ORB04950 
0RB04960 
ORB04970 
0RB04980 
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SUBROUTINE  RSWI JK( LAN , AL , INC , R , S , W , I , J , K) 

THIS  SUBROUTINE  CHANGES  FROM  RSW  COORDINATES  TO  UK  COORDINATES 


DOUBLE 
+   R31 
Rll  = 
R21  = 
R31  = 
R12  = 
R22  = 
R32  = 
R13  = 
R23  = 
R33  = 
I  =  Rl 
J  =  R2 
K  =  R3 
RETURN 
END 


PRECISION   INC,R,S,W,I,J,K,R11,R12,R13,R21,R22,R23, 
,R32,R33,LAN,AL 

DCOS(LAN)*DCOS(AL)  -  DSIN(LAN)*DCOS( INC)*DSIN( AL) 
DSIN(LAN)*DCOS(AL)  +  DSIN( AL)*DCOS(LAN)*DCOS( INC) 
DSIN(INC)*DSIN(AL) 

-DCOS(LAN)*DSIN(AL)-DSIN(LAN)*DCOS(INC)*DCOS(AL) 
-DSIN(LAN)*DSIN(AL)  +  DCOS(LAN)*DCOS( INC)*DCOS(AL) 
DSIN(INC)*DCOS(AL) 
DSIN(LAN)*DSIN(INC) 
-DCOS(LAN)*DSIN(INC) 
DCOS(INC) 

1*R  +  R12*S  +  R13*W 
1*R  +  R22*S  +  R23*W 
1*R  +  R32*S   +   R33*W 


SUBROUTINE  PQWRSW(TA , P ,Q , W ,R , S ,WN) 

*  THIS  SUBROUTINE  CHANGES  FROM  PQW  COORDINATES  TO  RSW  COORDINATES 

DOUBLE  PRECISION  P,Q,W,R, S ,WN,R11 ,R12,R13 ,R21 ,R22 ,R23 , 
+   R31,R32,R33,TA 
Rll  =  DCOS(TA) 
R12  =  DSIN(TA) 
R13  =  0.0 
R21  =  -DSIN(TA) 
R22  =  DCOS(TA) 
R23  =  0.  0 
R31  =  0.  0 
R32  =  0.0 
R33  =  1.0 

R  =  R11*P  +  R12*Q  +  R13---W 
S  =  R21*P  +  R22--'-Q  +  R23-'"W 
WN  =  R31*P  +  R32-Q  +R33-W 
RETURN 
END 

SUBROUTINE  RSWPQW( TA , R , S , W , P , Q , WN ) 

*  THIS  SUBROUTINE  CHANGES  FROM  RSW  COORDINATES  TO  PQW  COORDINATES 

DOUBLE  PRECISION  R,S ,W,P,QJWN,R11,R12 ,R13,R21 ,R22,R23, 
+   R31,R32,R33,TA 

Rll  =  DCOS(TA) 
R21  =  DSIN(TA) 
R31  =  0.  0 
R12  =  -DSIN(TA) 
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R22  =  DCOS(TA) 

R32  =  0.0 

R13  =  0.  0 

R23  =  0.  0 

R33  =  1.  0 

P   =  R11*R  +  R12*S   +  R13*W 

Q  =  R21*R  +  R22*S   +  R23*W 

WN  =  R31*R  +  R32*S  +R33*W 

RETURN 

END 

yr-iWrVrVcVrVr,-/r,)Vyryr-:V,;V->Vyf-sVyoV-;V'yr^^ 

STORE  ELEMENTS  IN  ARRAYS 

■fr^V-jV^V-;ryryrV«V-)V-!VycyryrV"V-^ 

SUBROUTINE  STORE (RI ,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY,TARAY,NUM, 
+    I)AP,AINRAY,APRAY,TT,TIMRAY) 

THIS  SUBROUTINE  STORES  THE  POSITION  AND  ELEMENTS  IN  ARRAYS  IN 
SINGLE  PRECISION  FORM,  FOR  PLOTTING 

DOUBLE  PRECISION  RI ,RJ,RK,R,TA, I , AP.TT 

DIMENSION  RIRAY(500) ,RJRAY(500) ,RKRAY(500) ,RARAY(500) ,TARAY(500) , 

+   AINRAY(500),APRAY(500),TIMRAY(500) 

RIRAY(NUM)  =  SNGL(RI) 
RJRAY(NUM)  =  SNGL(RJ) 

RKRAY(NUM)  =  SNGL(RK) 
RARAY(NUM)  =  SNGL(R) 
TARAY(NUM)  =  SNGL(TA) 
AINRAY(NUM)  =  SNGL(I) 
APRAY(NUM)  =  SNGL(AP) 
TIMRAY(NUM)  =  SNGL(TT) 
RETURN 
END 

*  INITIAL       POSITION,  VELOCITY 


SUBROUTINE  INPUTS ( RI , R J , RK , R , VI , VJ , VK , V , MU , QUIT , PI ) 

THIS  SUBROUTINE  GIVES  THE  USER  A  CHOICE  TO  EITHER  ENTER  THE 

INITIAL  POSITION  AND  VELOCITY  VECTOR  OR  TO  LET  THE  PROGRAM 

CALCULATE  THE  INITIAL  POSITION  AND  VELOCITY  FROM  USER  PROMPTED 

INPUTS 


SUBROUTINES  CALLED  FROM  THIS  SUBROUTINE: 
INELTS  =  Prompts  USER  FOR  ORBITAL  ELEMENTS 
IPOS   =  PROMPTS  USER  FOR  INITIAL  POSITION  (UK) 
IVEL   =  PROMPTS  USER  FOR  INITIAL  Velocity  (UK) 

DOUBLE  PRECISION  RI ,RJ,RK,R, VI , VJ,VK, V.MU.PI 
CHARACTERS,  QUIT 

PROMPT  USER  FOR  METHOD  TO  ENTER  INPUTS 
195   PRINT--'-',  '  IN"  WHICH  MANNER  WOULD  YOU  LIKE  TO  INPUT  THE  INITIAL* 
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POSITION  AND  VELOCITY  OF  THE  SATELLITE? 


1:  BY  Inputting  THE  INITIAL  POSITION  AND  VELOCITY' 
VECTORS  IN  THE  PERIFOCAL  COORDINATE  SYSTEM  fIJK)' 

2:  BY  LETTING  THE  SATELLITE  BE  PLACED  ON  THE  "I"1 
AXIS  OF  THE  (UK)  SYSTEM  AT  A  DESIRED  RADIUS  OF' 
PERIGEE(RP)  AND  INPUTTING  EITHER  A  DESIRED  RADIUS 
OF  APOGEE(RA),  A  DESIRED  ECCENTRICITY(E) ,  OR  THE' 


3: 

ENTER  1 


DESIRED  VELOCITY 
INCLINATION( I).  ' 
QUIT' 
2  OR  3:  ' 


AT  THAT  RADIUS,  AND  A  DESIRED 


PRINT* 

PRINT- 
PRINT* 
PRINT- 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
PRINT* 
READ*,ICHC 
PRINT*, ICHC 
CALL  EXCMS('CLRSCRN') 


USER  INPUTS  POSITION  AND  VELOCITY  VECTORS 
IF  (ICHC  . EQ.  1)  THEN 

CALL  IPOS(RI,RJ,RK,R) 

CALL  IVEL(VI,VJ,VK,V,R,MU) 

USER  INPUTS  ORBITAL  ELEMENTS  TO  GET  POSITION  AND  VELOCITY 
ELSEIF  (ICHC  .EQ.  2)  THEN 

CALL  INELTS(RI,RJ,RK,R,VI,VJ,VK,V,MU,PI) 

STOP  PROGRAM 

ELSEIF  (ICHC  .EQ.  3)  THEN 

QUIT  =  'N' 
ELSE 

PRINT*, 'INVALID  ENTRY!   TRY  AGAIN!' 

GOTO  195 
END  IF 
RETURN 
END 


.-  -.'.-  5>V  V.-  Vr  Vr  -,'.-  ■aftr  iV  ■ V  i;  V.-  Vr  i't  -.'r  •;'.-  Vr  -,V  -.V  -,'r  Vc  •>'?  Vr  Vr  -.',-  Vc  •>';  iV  •>'.-  1't  Vr  •>!  -.V  -,'r  -. 


SUBROUTINE  IPOS(RI ,RJ ,RK,R) 

THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  INITIAL  POSITION  OF  THE 

Satellite  IN  GEOCENTRIC-EQUATORIAL  COORDINATE  SYSTEM 

DOUBLE  PRECISION  RI,RJ,RK,R 

CHARACTER* 1,  CHOICE 
LOGICAL  CORREC 
CORREC   =  .FALSE. 

PROMPT  USER  FOR  VELOCITY  VECTOR 

180   IF(.  NOT. CORREC)  THEN 

CALL  EXCMS('CLRSCRN') 

PRINT*, 'ENTER  RADIUS  VECTOR  VALUES  IN  "KM"' 

PRINT*, 'RADIUS  OF  THE  EARTH  =  6400  KM' 

CORREC  =  .TRUE. 

PRINT*, 'ENTER  RI  : ' 

RE AD*, R I 

PRINT*, 'RI  =  ' ,RI,'KM' 

PRINT*, 'ENTER  RJ  :  ' 
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RE AD*, R J 

PRINT*, 'RJ  =  ' ,RJ,'KM' 

PRINT""-, 'ENTER  RK  :  ' 

READ*,RK 

PRINT-', 'RK  =  '  ,RK,'KM' 

CALCULATE  TOTAL  R 

R  =  DSQRT((RI**2)  +  (RJ**2)  +  (RK**2)) 

PRINT*, 'R  =  ' ,R, 'KM' 

IF  (R  .  LE.  6400.0)  THEN 

PRINT*, 'RADIUS  TO  SMALL! !   ENTER  NEW  VALUES! ! ' 

GOTO  180 
END  IF 

*        CHECK  WITH  USER  THAT  Values  ARE  CORRECT 

PRINT*,' ARE  THESE  VALUES  CORRECT?' 

PRINT'S1     ENTER  "Y"  OR  "N"  :' 

RE AD*, CHOICE 

CHOICE  =  'Y' 

PRINT*, CHOICE 

IF  (CHOICE. EQ. 'Y')  THEN 
CORREC   =  .TRUE. 

ENDIF 

GOTO  180 
ENDIF 
RETURN 
END 

SUBROUTINE  IVEL( VI , VJ , VK , V , R , MU) 

THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  INITIAL  VELOCITY  OF  THE 

Satellite 

DOUBLE  PRECISION  VI , VJ , VK, V,R, VCIR, VMAX.MU 

CHARACTER* 1,  CHOICE 
LOGICAL  CORREC 
CORREC  =  .FALSE. 

CALCULATE  ESCAPE  VELOCITY  AND  CIRCULAR  VELOCITY  AND  PROMPT  USER 
FOR  VELOCITY  VECTOR 
190   IF(.NOT. CORREC)  THEN 

CALL  EXCMS( 'CLRSCRN') 

VCIR  =  DSQRT(MU/R) 

VMAX  =  DSQRT((2.0*MU)/R) 

PR  I NT-S' CIRCULAR  VELOCITY  =  '  ,  VCIR,  '  KM/SEC' 

PR I NT*,* MAX I MUM  VELOCITY  =  ', VMAX, ' KM/ SEC' 

CORREC   =  .TRUE. 

PRINT*,' ENTER  VELOCITY  VECTOR  IN  (KM/SEC)' 

PRINT*, 'ENTER  VI  : ' 

RE AD*, VI 

PRINT*, 'VI  =  ', VI,' KM/SEC' 

PRINT*, 'ENTER  VJ  : ' 

READ---" ,  V J 
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PRINT*, 'VJ  =  ' ,VJ, 'KM/SEC' 

PRINT-, 'ENTER  VK  : ' 

READ*,VK 

PR I  NT* , ' VK  =  ' , VK , ' KM/ SEC ' 

CALCULATE  TOTAL  VELOCITY  (V) 

V  =  DSQRT((VI**2)  +  (VJ**2)  +  (VK**2)) 

PRINT*, 'V  =  ' ,V,'KM/SEC' 

CHECK  WITH  USER  THAT  VALUES  ARE  CORRECTS 

PR INT*,' ARE  THESE  VALUES  CORRECT?' 

PRINT*,1     ENTER  "Y"  OR  "N"  :  ' 

RE AD*, CHOICE 

CHOICE  =  'Y' 

PRINT*, CHOICE 

IF  (CHOICE. EQ.  *Yf )  THEN 

CORREC   =  .TRUE. 
END  IF 
IF  (V  . GE.  VMAX)  THEN 

PRINT*,' VELOCITY  IS  GREATER  THAN  THE  ESCAPE  VELOCITY!! ' 

PRINT*, 'RE-ENTER  VELOCITY!  !  !  ' 

CORREC  =  .FALSE. 
END  IF 
GOTO  190 
END  IF 
RETURN 
END 

SUBROUTINE  INELTS(RI,  RJ,  RK,  R,  VI,  VJ,  VK,  V,MU,PI) 

SATELLITE  PLACED  ON  ' I '  AXIS  AND  USER  SUPPLY  ORBITAL  ELEMENTS  TO 

GET  INITIAL  POSITION  AND  VELOCITY 

DOUBLE  PREC I S ION  RI , RJ , RK , R , VI , V J , VK , V , MU , I , ENR , A , E , RP , RA , PI , VMAX 
CHARACTER* 1, CHOICE 

PROMPT  USER  FOR  PERIGEE  RADIUS 

PRINT*, 'ENTER  RADIUS  OF  PERIGEE(RP)  IN  (KM),  FOR  EXAMPLE:' 

PRINT*, 'LOW  EARTH  ORBIT  (LEO),  RP  =  6600.0  KM' 

PRINT*, 'GEOSYNCROUNOUS  ORBIT,  RP  =  42241.1  KM' 

PRINT*, 'ENTER  RP:  ' 

PRINT'S '"RP"  MUST  BE  >  6400KM' 

READ*,RP 

PRINT*, RP 

CHECK  FOR  VALID  RADIUS 
IF  (RP  . LT.  6400. 0)  THEN 

PRINT*,'YOUR  "RP"  IS  TO  SMALL!! ' 

GOTO  198 

END  IF 

PROMPT  USER  FOR  TYPE  OF  INPUT 

PRINT*, 'DO  YOU  WANT  TO  ENTER  THE  ECCENTRICITY  (E),  ' 
PRINT*, 'RADIUS  OF  APOGEE  (RA),  OR  VELOCITY  (V)?' 
PRINT*, 'ENTER  "E",  "R" ,  OR  "V":  ' 
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READ*, CHOICE 
PRINT''-,  CHOICE 

CALL  EXCMS('CLRSCRN') 

USER  ENTERS  Eccentricity  AND  SEMI -MAJOR  AXIS,  ENERGY  AND  VELOCITY 
IS  CALCULATED  IN  THAT  ORDER 
IF  (CHOICE  . EQ.  'E' )  THEN 

PRINT*, 'ENTER  ECCENTRICITY  (E):' 

PRINTS, '0.  0   <=  E   <    1.0' 

READ*,E 

PRINT'S  E 

CHECK  FOR  VALID  ECCENTRICITY 

IF  ((E  .  LT.  0.0)  .OR.  (E  .  GE.  1.0))  THEN 

PRINT'S  '  INVALID  "E"1 

GOTO  198 
END  IF 

A  =  RP/(1-E) 
ENR  =  -ML'/ (2.  0*A) 

V  =  DSQRT(2*(ENR+(MU/RP))) 

USER  INPUTS  RADIUS  OF  APOGEE  AND  ECCENTRICITY  IS  CALCULATED 
THEN  SEMI-MAJOR  AXIS,  ENERGY  AND  THEN  VELOCITY. 
ELSEIF  (CHOICE  .  EQ.  'R')  THEN 

PRINT'S  ' ENTER  RADIUS  OF  APOGEE  (RA)  IN  KM:  ' 
PRINT'S '"RA"  MUST  BE  >="RP",  "RP"  =  '  ,RP 
READ* , RA 
PRINTS  RA 

CHECK  FOR  VALID  RADIUS  OF  APOGEE 
IF  (RA  .LT.  RP)  THEN 

PRINT'S 'YOUR  "RA"  IS  TO  SMALL!!' 

GOTO  198 
END  IF 

E  =  (RA-RP)/(RA+RP) 
A  =  RP/O-E) 
ENR  =  -ML'/ (2.  0*A) 

V  =  DSQRT(2*(ENR+(MU/RP))) 

USER  INPUTS  MAGNITUDE  OF  VELOCITY,  PROGRAM  PROVIDES  CIRCULAR 
AND  ESCAPE  VELOCITY  FOR  COMPARISON  AND  TO  CHECK  FOR  VALID 
INPUTS 
ELSEIF  (CHOICE  . EQ.  'V')  THEN 

PRINT'S 'ENTER  VELOCITY  IN  KM/SEC:  ' 

PRINT'S 'THE  MINIMUM  VELOCITY  ALLOWED  IS  FOR  A  CIRCULAR  ORBIT* 

VCIRC  =  SQRT(SNGL(MU/RP)) 

PRINT'S' ORBIT.  V(  Circular)  =  *  ,  VCIRC,'  KM/S' 

VMAX  =  DSQRT(2*(MU/RP)) 

PR  I  NTS' THE  MAXIMUM  VELOCITY  <  ',VMAX,'  KM/S' 

READ*,V 

PRINT'S  V 

IF  (V  .LT.  VCIRC)  THEN 

PRINT*, 'VELOCITY  TO  SMALL!' 

GOTO  198 
END  IF 
IF  (V  . GE.  VMAX)  THEN 
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PRINT*, 'VELOCITY  TO  GREAT! ! ' 
GOTO  198 
ENDIF 

ELSE 

PRINT*, 'INVALID  ENTRY!  TRY  AGAIN' 

GOTO  198 
ENDIF 

*  INCLINATION  NEEDED  TO  GIVE  Velocity  A  Direction 
PRINT*, 'ENTER  INCLINATION  (I)  IN  DEGREES:' 
READ* , I 

PRINT*, I 

I  =  (PI/180. 0)*I 
VK  =  V*DSIN(I) 
VJ  =  V*DCOS(I) 
VI  =  0.0 

*  RADIUS  VECTOR  SET 
RI  =  RP 

RJ  =  0.  0 
RK  =  0.  0 
R  =  RP 
RETURN 
END 

*  CALCULATE  THE  ORBITAL  ELEMENTS 


SUBROUTINE  CALCEL(RI ,RJ,RK,R, VI , VJ, VK,V,EI ,EJ,EK,E , A, I ,LAN, 
+  LP,TA,PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

THIS  SUBROUTINE  CALLS  THE  INDIVIDUAL  SUBROUTINES  TO  CALCULATE  THE 
ORBITAL  ELEMENTS 

THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES( RETURNED  VALUES) 

ENERGY  =  ENERGY  PER  MASS  (ENR) 

ANGMOM  =  ANGULAR  MOMENTUM  (H,HI ,HJ,HK) 

NODE   =  NODE  VECTOR  (N,NI,NJ,NK) 

LATREC  =  SEMI-LATUS  RECTUS  (P) 

ECC    =  ECCENTRICITY'  (E,EI,EJ,EK) 

SMAXIS  =  SEMI -MAJOR  AXIS  (A) 

INCL   =  INCLINATION  (I) 

ASNODE  =  LONGITUDE  OF  ASCENDING  NODE  (LAN) 

ARP    =  ARGUMENT  OF  PERIGEE  (AP) 

IJKPQW  =  'UK'  SYSTEM  TO  'PQW'  SYSTEM 

TANOM  =  TRUE  ANOMALY  (TA) 

ARLAT  =  ARGUMENT  OF  LATITUDE  (AL) 

LONPER  =  LONGITUDE  OF  Perigee  (L?) 

TLON   =  TRUE  LONGITUDE  (TL) 

PERIOD  =  PERIOD  (PER) 

ECC AN  =  ECCENTRIC  ANOMALY  (EA) 

MEANMO  =  MEAN  MOTION  (MM) 

MEANAN  =  MEAN  ANOMALY  (MA) 

TFLGKT  =  TIME  OF  FLIGHT  (TF) 

DOUBLE  PRECISION  RI ,RJ,RK,R, VI , VJ,VK, V,EI ,EJ,EK,E , A, I ,LAN, AL, 
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ORB08360 
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ORB08430 
ORB08440 
ORB08450 
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ORB08470 
0RB08480 
0RB08490 
ORB08500 
ORB08510 
ORB08520 
ORB08530 
ORB08540 
ORB08550 
ORB08560 
ORB08570 
ORB08580 
ORB08590 
ORB08600 
ORBO8610 
ORB08620 
ORB08630 
ORB08640 
ORB08650 
ORB08660 
ORB08670 
ORB08680 
ORB08690 
ORB08700 
ORB08710 
ORB08720 
ORB08730 
ORB08740 
ORB08750 
ORB08760 
ORB08770 
ORB08780 
ORB08790 
ORB08800 
0RB08810 
ORB08820 
ORB08830 
ORB08840 
ORB08850 
ORB08860 
ORB08870 
ORB08880 
ORB08890 
ORB08900 
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+   LP.TA.PER.EA.MA.AP.TF.HI.HJ.HK.H.NI.NJ.NK.N.P.PI.MU.MM.ENR, 
+   TL,RP,RQ,RW,NP,NQ,NW 

CALL  ENERGY(V,R,MU,ENR) 

CALL  ANGMOM( RI , RJ , RK , VI , VJ , VK , HI , HJ , HK , H) 

CALL  NODE(HI,HJ,NI,NJ,NK,N) 

CALL  LATREC(H,P,MU) 

CALL  ECC(RI,RJ,RK,R,VI,VJ>VK,V,EI>EJ,EK)E>MU) 

CALL  SMAXIS(MU,ENR,A) 

CALL  INCL(KK,H,I,PI) 

SPECIAL  CASE  IF  INCLINATION  =0.0 
IF  (I.NE.O.  0)  THEN 

CALL  ASNODE(NI,N\LAN,NJ,PI) 

CALL  ARP(NI,NJ,N,EI,EJ,EK,E,AP,PI,NP,NQ,LAN) 

ELSE 

LAN  =  0.  0 

A?  =  0.  0 
ENDIF 

COORDINATE  TRANSFORMATION  OF  'R1  AND  'V'  VECTORS 

CALL  IJKPQW(LAN,AP,I,RI,RJ,RK,RP,RQ,RW) 

CALL  IJKPQW(LAN,AP,I,NI,NJ,NK,NP,NQ,NW) 

CALL  TANOM( E I , E J , EK , E , RI , R J , RK , RP , RQ , RW , R , V I , V J , VK , TA , P I ) 

SPECIAL  CASE  FOR  Inclination  =0.0 
IF  (I  . NE.  0. 0)  THEN 

CALL  ARLAT( NI , N J , NK , N , RI , R J , RK , R , AL , P I , TA , AP ) 
ELSE 

AL  =  TA 
ENDIF 

CALL  LONPER(LAN,AP,LP) 
CALL  TLON(LAN,AP,TA,TL) 
CALL  PERIOD(A,PER,PI,MU) 
CALL  ECCAN(E,TA,EA,PI) 
CALL  MEANMO(A,MM,MU) 
CALL  MEANAN(EA,E,MA) 
CALL  TFLGHT(MM,MA,TF) 
RETURN 
END 


.'.  JL.  »*-  -'-  -'-  J—J-  J.  J. 
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SUBROUTINE  ENERGY( V , R , MU , ENR) 

THIS  SUBROUTINE  CALCULATES  THE  ENERGY  OF  THE  ORBIT 

DOUBLE  PRECISION  V,R,MU,ENR 

ENR  =  ((V**2)/2)  -  (MU/R) 

RETURN 

END 

SUBROUTINE  ANGMOM(RI , RJ , RK , VI , VJ , VK , HI , HJ , HK , H) 
THIS  SUBROUTINE  CALCULATES  THE  ANGULAR  MOMENTUM 
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DOUBLE  PRECISION  RI ,RJ,RK, VI , VJ, VK,HI ,HJ ,HK,H 

KI  =  (RJ  *  VK)  -  (RK  *  VJ) 
HJ  =  (RK  *  VI)  -  (RI  ••••'  VK) 
HK  =  (RI  *  VJ)  -  (RJ  *  VI) 
H  =  DSQRT((HI**2)  +  (HJ**2)  +  (HK**2)) 

RETURN 
END 

SUBROUTINE  NODE(HI ,HJ,NI ,NJ,NK,N) 
*     THIS  SUBROUTINE  CALCULATES  THE  NODE  VECTOR 

DOUBLE  PRECISION  HI ,HJ,NI ,NJ,NK,N 

NI  =  -HJ 

NJ  =  HI 

NK  =  0.0 

N  =  DSQRT((NI**2)  +  (NJ**2)) 

RETURN 
END 

SUBROUTINE  LATREC(H,P ,MU) 

THIS  SUBROUTINE  CALCULATES  THE  SEMI-LATUS  RECTUM 

DOUBLE  PRECISION  H,P,MU 

P  =   (H**2)/MU 

RETURN 
END 


J.  «'-  Jm  -' .  Jm  Jm  J<  , ' .  - '.  JU  »'.  J.  -•  -  *3m  ~.  Jm  - '-  Jm  - '  -  J -  •'-  J  -  J -  -J.  -'  - 
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SUBROUTINE  ECC(RI ,RJ,RK,R, VI , VJ, VK, V,EI ,EJ,EK,E ,MU) 
THIS  SUBROUTINE  CALCULATES  THE  ECCENTRICITY 

DOUBLE  PRECISION  RI ,RJ,RK,R, VI , VJ, VK, V,EI ,EJ ,EK,E ,MU,D0T 


A 


CALCULATE  DOT  PRODUCT  OF  'R'  AND  'V'  VECTORS 

DOT  =  (RI*VI)  +  (RJ*VJ)  +  (RK*VK) 

EI  =  (1.0D+00/MU)  *  (((V**2)  -  (MU/R))  *  RI  -  (DOT)*VI) 

EJ  =  (1.0D+00/MU)  *  (((V**2)  -  (MU/R))  *  RJ  -  (DOT)-VJ) 

EK  =  (1.0D+00/MU)  *  C((V**2)  -  (MU/R))  *  RK  -  (DOT)*VK) 

E  =  DSQRT((EI**2)  +  (EJ**2)  +  (EK**2)) 

RETURN 

END 


SUBROUTI NE  SMAXI S ( MU , ENR , A ) 

THIS  SUBROUTINE  Calculates  THE  SEMI -MAJOR  AXIS 
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DOUBLE  PRECISION  MU,ENR,A 

A  =  -MU/(2*ENR) 
RETURN 

END 

SUBROUTINE  INCL(HK,H , I ,PI) 
*     THIS  SUBROUTINE  CALCULATES  THE  INCLINATION 
'I'  ALWAYS  LESS  THAN  180  DEGREES 

DOUBLE  PRECISION  HK,H,I,PI 

I  =  DACOS(HK/H) 

RETURN 
END 

VrVriVVfV.-VrVrVr^VVrVrVrVr'jV-j'.-Vr'j't'Vr-V-.VV.-VcVrVrVryr'j'f-A-'s'rVf  ArriVVcVrVrVrVriV^VryfiVyrVrVrVriVVrVr,)ViV-iV-,VVf')V'5V';V'A,V?';V-.'f-,V--'.-'!V 

SUBROUTINE  ASNODE ( NI , N , LAN , NJ , PI ) 

THIS  SUBROUTINE  CALCULATES  THE  LONGITUDE  OF  THE  ASCENDING  NODE 

IF  'NJ'  >  0       THEN  'LAN'  <  180  DEGREES 

DOUBLE  PRECISION  NI,N,LAN,NJ,PI 

LAN  =  DATAN2(NJ,NI) 
IF  (LAN  .  LT.  0.  0)  THEN 

LAN  =  (2---PI)  +  LAN 
END  IF 
RETURN 
END 
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SUBROUTINE  ARP(NI ,NJ ,N ,EI ,EJ ,EK ,E , AP ,PI ,NP ,NQ,LAN) 
THIS  SUBROUTINE  CALCULATES  THE  ARGUMENT  OF  Perigee 
IF  'EX'  GREATER  THAN  0  THEN  'AP'  <  180 
VARIABLE  TEMP  USED  AS  A  Temporary  VALUE  FOR  ARCTAN 

DOUBLE  PRECISION  NT  ,NJ,N,EI ,EJ,EK,E , AP ,PI ,NQ, NT, TEMP, LAN 

IF  ((EI  .  EQ.  0.0)  .AND.  (EJ  .  EQ.  0.0))  THEN 

AP  =  0.  0 
ELSE 

TEMP  =  DATAN2(EJ,EI) 
IF  (TEMP  . GT.  LAN)  THEN 

AP  =  TEMP  -  LAN 
ELSE 

AP  =  (2*PI)  -  (LAN  -  TEMP) 
END  IF 
IF  (  AP  .LT.  0.0)  THEN 

AP  =  (2*PI)  +  AP 
END  IF 

IF  (AP  .GT.  (2*PI))  THEN 
AP  =  A?  -  (2*PI) 
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ORB10560 
ORB10570 
ORB10580 
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DIF 

::- 

RET  - 


-■»  J(a  y-  JL  -'»  JU  Jm  Jg 


SUBROUTINE  TANGM( EI ,E J,EK,E ,RI ,RJ,RK,RP ,RQ, RW ,R, VI , VJ , VK , 

+  TA.PI) 

THIS  SUBROUTINE  CALCULATES  THE  TRUE  Anomaly 
IF  (R  DOT  V)  >  0  THEN"  TA  <  130  DEGREES 

DOUBLE  PRECISION  DOT,EI,EJ,EK,E,RI,RJ,RK,R,VI,VJ,VK,TA,PIJ 

+   rp  ,rq  ,:-.;; 

TA  =  DATAN2(RQ,RP) 
IF  (TA  .LT.  0.  0  )  THEN 
TA  =  (2  *  PI)  +  TA 
END  IF 
RETURN 
END 


-»'-  ->«•.«*..■-..'.  -■-  JL 


J—X^Xm^m^mi 


SUBROUTINE  ARLAT(  NI ,  NJ ,  NK  ,N .  ?  I  .  ?  J  .  ?  E  .  ?. ,  AL,PI  ,TA,  A? 
THIS  SUBROUTINE  CALCULATES  TEE  ARGUMENT  IE  LATITUDE 

IF  (RK  >  0)  TEEN  AL  <  15:  EECEEES 

DOUBLE  3EECISI0N  NI,NJ,NK,K,RI,RJ,BK,R,AL,PI,TA,AI 


RETURN 

END 





SUBROUTINE  LONPER(LAN, AP ,LP) 

THIS  SUBROUTINE  CALCULATES  THE  LONGITUDE  3F  PERIGEE 

DOUBLE  PRECISION  LAN,AP,LP 

LP  =  LAN  +  A? 
RETURN 

END 


ikkkkXXXXXkkkki 


SUBROUTINE  TLON( LAN , AP , TA . TE 

THIS  SUBROUTINE  CALCULATES  TEEE  TRUE  LONGITUDE  AT  ZESCE 

DOUBLE  PRECISION  LAN,AP,TA,TL 

TL  =  A?  +  LAN  +  TA 
RETURN 

END 


2  10590 
3RB10600 

31061C 
)RB1062C 
)RB10630 
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DRB10650 

ORB10660 
DRB10670 
". -310680 
3RB10690 
..-310700 
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ORB10740 
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DRB1079C 
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3RB1C  31C 
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SUBROUTINE  PERIOD( A , PER , PI , MU) 

*  THIS  SUBROUTINE  CALCULATES  THE  PERIOD 

DOUBLE  PRECISION  A, PER, PI, MU 

PER  =  2.0D+00*(PI)*DSQRT((A**3)/MU) 

RETURN 
END 

SUBROUTINE  ECCAN(E ,TA,EA,PI) 

*  THIS  SUBROUTINE  CALCULATES  THE  ECCENTRIC  Anomaly 

DOUBLE  PRECISION  E,TA,EA,PI 

EA  =  DACOS((E  +  DCOS(TA))/(l.  OD+00  +  E*DCOS(TA) ) ) 
IF  (TA  .GT.  PI)  THEN 
EA  =  (2*PI)  -  EA 

END  IF 

RETURN 

END 

SUBROUTINE  MEANMO( A, MM ,MU) 

THIS  SUBROUTINE  CALCULATES  THE  MEAN  MOTION 

DOUBLE  PRECISION  A,MM,MU 

MM  =  DSQRT(MU/(A**3)) 

RETURN 

END 
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SUBROUTINE  MEANAN(EA ,E ,MA) 

THIS  SUBROUTINE  CALCULATES  THE  MEAN  Anomaly 

DOUBLE  PRECISION  EA,E,MA 

MA  =  EA  -  E*DSIN(EA) 

RETURN 

END 

yoVyr^yrycVcyrVfyryoV*}VyoWoVyr^VrVryriWr*Tfcyr*:fr^^^ 

SUBROUTINE  TFLGHT(MM,MA,TF) 

*  THIS  SUBROUTINE  CALCULATES  THE  TIME  OF  FLIGHT 

DOUBLE  PRECISION  MM,MA,TF 
TF  =  (1/MM)*MA 
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RETl 
E 








it 


CALCULATE  UNPERTURBEE  :~BI7 


SUBROUTINE  UNPRET( DT , PER , AL , LAN , AP , I , RI , RJ , RK , R , 
+   VI ,VJ,VK,V ,MU, PI,H,A,E3 N ,TA ,P  MM  MA  EA 

+   77  \  7 ,  NUM ,RIRAY  ,  77ZAY  ]  Z77AY  .  ZA7AY  .  7A7AY  ,  AIN7AY  ,  A?7A7  7IMRAY  . 
+   77 

THIS  SUBROUTINE  CALCULATE  THE  UNPERTURBEE  DRBIT 

THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 
NEWELT  =  CALCULATE  NEW  ELEMENTS  A77ER  TIME  S7Z? 
NEWPOS  =  CALCULATE  NEW  POSITION  A777R  TIME  STEP 
NEWVEL  =  CALCULATE  NEW  VELOCITY  A77ER  TIME  STEP 

S7CPE  =  stores  position  in  a?pays 

double  precision  t,dt,per,al,lan,aj . [,ri,rj,rk,r,vi,vj,vk,vj 
+  mu,pi,h,a,e,n,ta,p,mm,ma,ea,tf,tt 

dimension  raray(500),taray(500),riray(500),rjray(500), 

+  RKRAY(500),AINRAY(500),APRAY(500),T::'.RAY  500) 


SET  TRUE  ANOMALY  TO  NEGATIVE  SC  LOOP 

IF  (TA  .  GT.  6.  21)  THEN 
7A  =  7A  -  (2*PI) 

END  17 


:an  be  ex 


to: 


CONTINUE  AROUND  ORB  17  71 LI 
230   IF  ((TA  .  LE.  6.  21)  .AND.  ( 


VSI"  7:  PERIGEE 
1Z. FEE'   THEN 


*  Increment  TRUE  TIME 
TT  =  TT  +  DT 

CALL  NEWELT(  MM  ,  MA  .  Z  .  ZA  ,  7A  .  77  .  17  .  ?  I  .  PER) 
CALL  NPOS(RI,RJ,RK,R,LAN,AP,I,  TA,A,E) 
CALL  NVEL(  E  ,  P  , TA ,  LAN ,  AF  ,  I ,  VI ,  VJ ,  VK  .  V  ,  .TO) 

*  INCREMENT  STEP  COUNTER  AND  S70RE  VALUES 
NU'M  =  MUM  +  1 

CALL  ST0RE(RIJRJJRK,R,TA,RIRAY,RJRAY,R7?V.V, 
+  RARAY , TARAY , NUM , I , AP , AINRAY , APRAY , 

+  TT,TIMRAY) 

*  INCREMENT  TIME  STEP  COUNTER 
T=  T  +  DT 

GOTO  230 
END  IF 
RETURN 
END 

*  CALCULATE  THE  UNPERTURBED  NEW  ELEMENTS 
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SUBROUTINE  N:EWELT(  MM ,  MA ,  E  ,  EA , TA , TF ,  DT ,  PI ,  PER ) 

THIS  SUBROUTINE  CALCULATES  THE  Unperturbed  NEW  ELEMENTS 

THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 
NEA  =  NEW  ECCENTRIC  ANOMALY 
NTA  =  NEW  TRUE  ANOMALY 

DOUBLE  PRECISION  MM,MA,E ,EA,TA,TF,DT,PI ,PER 

*  Increment  TIME  OF  FLIGHT  AND  CHECK  IF  TF  GREATER  THAN  PERIOD 

TF  =  TF  +  DT 

IF  (TF  .GT.  PER)  THEN 

TF  =  TF  -  PER 
END  IF 

CALCULATE  MEAN  ANOMALY  AND  USE  TO  FIND  ECCENTRIC  Anomaly  THEN  NEW 

*  TRUE  ANOMALY 
MA  =  MM*(TF) 
CALL  NEA(MA,E,EA) 
CALL  NTA(EA,E,TA,PI) 
RETURN 

END 

*  CALCULATE  PERTURBED  ORBIT 

SUBROUTINE  PRETUR( DT , PER , AL , LAN , AP , I , RI , R J , RK , R , 
+   VI.VJ.VK.V.FR.FS.FW.MU.P^H.A.E^.TA.P.MM.MA.EA, 
+   TF , T , NUM , RIRAY , R JRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY , TIMRAY , 
+   TT,TFEA,TFSU,TFMO,TFDRA,TDI,TDA,TDE,TDMM,TDMA,TDLAN,TDH,TDAP) 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBED  ORBIT. 

THIS  SUBROUTINE  CALLS  THE  FOLLOWING  SUBROUTINES: 

TFORCE  =  CALCULATE  THE  TOTAL  PERTURBING  FORCE  ON  THE  SATELLITE 

PNEWEL  =  CALCULATE  THE  Perturbed  NEW  ELEMENTS 

NPOS   =  NEW  POSITION  AFTER  TIME  STEP 

*  NVEL   =  NEW  VELOCITY  AFTER  TIME  STEP 
PERIOD  =  PERIOD  OF  PERTURBED  ORBIT 

STORE   =  STORE  POSITION  AND  ELEMENTS  IN  ARRAYS  FOR  PLOTTING 

DOUBLE  PRECISION  T,DT,PER, AL,LAN, AP, I ,RI ,RJ,RK,R ,VI ,VJ, VK,V, 
+  FR,FS,FW,MU,PI,H,A,E,N,TA,P,MM,MA,EA,TF,TT, 

+   DI,DA,DE,DMM,DMA,DLAN,DH,DAP,EI,EJ,EK,HI,HJ,LP,M, 
+   DVR,DVS,DVW,DVI,DVJ,DVK 

DIMENSION  RARAY(500) ,TARAY(500) ,RIRAY(500) ,RJRAY(500) , 
+  RKRAY(500),AINRAY(500),APRAY(500),TIMRAY(500) 

*  SET  MEAN  RADIUS  OF  EARTH 
RE  =  6400.  0 

DT  =  PER/50 
T  =  DT 

IF  (TA  .GT.  6.  21)  THEN 
TA  =  TA  -  (2*PI) 
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END  IF 

IF  (TF  .  GE.  PER)  THEN 

TF  =  TF  -  PER 
END  IF 

CONTINUE  Around  ORBIT  FOR  ONE  PERIOD 

240  IF  ((TF  .  LT.  PER)  .AND.  (T  .  LT.  PER))  THEN 

•        INCREMENT  TRUE  TIME 
TT  =  TT  +  DT 

CALL  TFORCE(AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+  TT,FR,FS,FW,MU,PI, 

+  FEA,FSU,FMO,FDRA,FOR, 

+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,TF,P, 

+  MM,N,H,HI,HJ,DT) 

CALL  PNEWEL(FR,FS,FW,H,R,A,E,N,TA,DT,I,LAN,AL, 
+  AP,P,MM,MA,EA,TF,T,MU,PI, 

+  DI,DA,DE,DMM,DMA,DLAN,DH,DAP) 

CALL  NPOS(RI,RJ,RK,R,LAN,AP,I,  TA,A,E) 
CALL  NVEL(EJP,TA,LANJAP,I,VI,VJ,VK)V,MU) 

CALCULATE  NEW  PERIOD  AND  RESET  TIME  STEP  AND  TIME  COUNTER 

■  IF  NOT  AT  END  OF  ORBIT 

IF  (T  .LT.  (PER-DT))  THEN 

CALL  PERIOD(A,PER,PI,MU) 

DT  =  PER/50 

T  =  TF 
END  IF 

■  INCREMENT  STEP  COUNTER 
NUM  =  NUM  +  1 

241  CALL  STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY, 

+  RARAY , TARAY , NUM , I , AP , AINRAY , APRAY , 

+  TT,TIMRAY) 


TOTAL 
TDI  = 
TDA  = 
TDE  = 
TDMM  = 
TDMA  = 
TDLAN 
TDH  = 
TDAP  = 
TFEA  = 
TFSU  = 
TFMO  = 
TFDRA 


ELEMENT  CHANGES 
TDI  +  SNGL(ABS(DI)) 
TDA  +  SNGL(ABS(DA)) 
TDE  +  SNGL(ABS(DE)) 
■■   TDMM  +  SNGL(ABS(DMM)) 
'  TDMA  +  SNGL(ABS(DMA)) 
=  TDLAN  +  SNGL(ABS(DLAN)) 


TDH  + 
■  TDAP 
=  TFEA 
=  TFSU 
=  TFMO 


=  TFDRA 


SNGL(ABS(DH)) 
+  SNGL(ABS(DAP)) 

+  FEA 
+  FSU 
+  FMO 

FDRA 


CHECK  FOR  IMPACT 
IF  (R  . LE.  RE)  THEN 

PRINT*, 'SATELLITE  WILL  IMPACT  THE  EARTH!! ' 

T  =  PER 
END  IF 

INCREMENT  TIME  COUNTER 
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T  =  T  +  DT 

GOTO  240 
END  IF 
RETURN 
END 

*  CALCULATE  THE  PERTURBING  FORCES 


SUBROUTINE  TF0RCE( AL,LAN, AP , I ,RI ,RJ,RK,R, VI , VJ, VK, V,TT, 
+  FR,FS,FW,MU,PI,FEA,FSU,FMO,FDRA,FOR, 

+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,TF,P, 

+  MM,N,H,HI,HJ,DT) 

THIS  SUBROUTINE  SUMS  ALL  THE  PERTURBING  FORCES  FOR  THE  TOTAL 

PERTURBING  FORCE. 


* 


Vc 


* 
* 


THE  FOLLOWING  SUBROUTINES  WERE  CALLED: 
OBERT  =  OBLATENESS  OF  THE  EARTH 
FSUN   =  GRAVITATIONAL  Attraction  OF  THE  SUN 
FMOON  =  GRAVITATIONAL  Attraction  OF  THE  MOON 
FDRAG  =  DRAG  FORCES 

DOUBLE  PRECISION  FER,FES ,FEW,FSR,FSS ,FSW,FMR,FMS ,FMW,MU,PI , 
+   FDR,FDS,FDW,FR,FS,FW,RI,RJ,RK,R,AL,I,TT,LAN,AP,VI,VJ,VK,V, 
+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,TF,P, 

+  MM,N,H,HI,HJ,DT 

CALL  OBEARTC RI , RJ , RK , R , AL , I , FER , FES , FEW , MU) 

CALL  FSUN(TT,RI,RJ,RK,R,FSR,FSS,FSW,PI) 

CALL  FMOON(TT,RI,RJ,RK,R,FMR,FMS,FMW,PI) 

CALL  FDRAG(RI,RJ,RK,R,VI,VJ,VK,V,LAN,AP,I,FDR,FDS,FDW, 
+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 

+  MM,N,H,HI,HJ,DT) 

SUM  VECTOR  FORCES 
FR  =  FER  +  FSR  +  FMR  +  FDR 
FS  =  FES  +  FSS  +  FMS  +  FDS 
FW  =  FEW  +  FSW  +  FMW  +  FDW 


CALCULATE  TOTAL  FORCE  FROM  EACH,  AND  TOTAL  OF  ALL 
FEA  =  SNGL(SQRT((FER**2)+(FES**2)+(FEW**2))) 
FSU  =  SNGL(SQRT((FSR**2)+(FSS**2)+(FSW** 2))) 
FMO  =  SNGL(SQRT((FMR**2)+(FMS**2)+(FMW**2))) 
FDRA  =  SNGL(SQRT((FDR**2)+(FDS**2)+(FDW**2))) 
FOR  =  SNGL(SQRT((FR**2)+(FS**2)+(FW**2))) 

RETURN 
END 

yryrVriWfy?V?VfV>-Vr^VVoVycyoVV??V-,yywwryr^ 

SUBROUTINE  OBEART(RI ,RJ,RK,R , AL, I , FER, FES ,FEW ,MU) 
*     THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  THE 
OBLIQUENESS  OF  THE  EARTH. 
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DOUBLE  PRECISION  J2 ,RE ,FER,FES ,FEW,RI ,RJ,RK,R, AL, I ,MU,M 

J2  =  1. 082364D-03 
RE  =  6. 3762D+03 

FER  =  ((-3. 0D+00*MU*J2*(RE**2))/(2.  0D+00*(R**4)))* 
+      (l.OD+OO  -  (3.0D+00*  ((DSIN(I))**2)  *  (  (DSIN(  AL)  )**2)  ) ) 
FES  =  ((-3. 0D+00*MU*J2*(RE**2))/(R**4))* 

+      (((DSIN(I))**2)*  (DSIN(AL))*(DCOS(AL))) 

FEW  =  ( ( -  3 . 0D+00*MU* J2*( RE**2 ) ) / ( R**4 ) ) * 
+      (DSIN(I)*DCOS(I)*DSIN(AL)) 

RETURN 
END 

SUBROUTINE  FSUN(TT,RI ,RJ,RK,R,FSR,FSS ,FSW,PI ) 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  THE  SUN 

THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

SUNPOS  =  SUNS  POSITION  ORBITING  AROUND  EARTH 

*  HEVBOD  =  PERTURBING  FORCE  FROM  A  Heavenly  BODY 

DOUBLE  PRECISION  FSR,FSS ,FSW, PI , 
+    RSI,RSJ,RSK,SLAN,SI,SAL,SMU,TT,RI,RJ,RK,R,RS 

SUNS  GRAVITATIONAL  PARAMETER 

SMU  =  1.  3271544D+11 

CALL  SUNPOS(TT,RSI,RSJ,RSK,RS,SLAN,SI,SAL,PI) 

CALL  HEVBOD(RI,RJ,RK,R,RSI,RSJ,RSK,RS,SLAN,SAL,SI,SMU,FSR,FSS,FSW 

RETURN 

END 

SUBROUTINE  FMOON(TT,RI ,RJ,RK,R,FMR,FMS ,FMW,PI) 

THIS  SUBROUTINE  CALCULATES  THE  PERTURBING  FORCE  DUE  TO  The  MOON 

*  THE  FOLLOWING  SUBROUTINE  ARE  CALLED: 

MONPOS  =  MOONS  POSITION  ORBITING  AROUND  THE  EARTH 
HEVBOD  =  PERTURBING  FORCE  FROM  A  HEAVENLY  BODY 

DOUBLE  PRECISION  FMR.FMS ,FMW,RMI ,RMJ,RMK,MLAN,MI ,MAL,MMU, 
+  TT,RI,RJ,RK,R,RM,PI 

*  MOONS  GRAVITATIONAL  PARAMETER 
MMU  =  4. 90287D+03 

CALL  MONPOS(TT,RMI,RMJ,RMK,RM,MLAN,MI,MAL,PI) 

CALL  HEVBOD( RI , RJ , RK , R , RMI , RMJ , RMK , RM , MLAN , MAL , MI , MMU , FMR , FMS , FMW 

RETURN 

END 

SUBROUTINE  HEVBOD( RI , RJ , RK , R , RPI , RP J , RPK , RP , LAN , AL , INC , MUP , 
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+  FHR.FHS.FHW) 

THIS    SUBROUTINE   CALCULATES   THE   PERTURBING   FORCE   DUE   TO   A 

HEAVENLY  BODY. 

THE   FOLLOWING   SUBROUTINE  WAS   CALLED: 
IJKRSW  =    'UK'    SYSTEM  TO  THE    'RSW'    SYSTEM 

DOUBLE   PRECISION     DOT,FHI ,FHJ,FHK,RI ,RJ,RK,R,RPI ,RPJ,RPK,RP, 
+        LAN'.AL.INC.MUP.I.J.K.IP.JP.KP.M.FHR.FHS.FHW 

CALCULATE  UNIT  VECTOR  FOR  SATELLITE  AND  PERTURBING  BODIES   POS 

I   =  RI/R 

J  =  RJ/R 

K  =  RK/R 

IP  =  RPI/RP 

JP  =  RPJ/RP 

KP  =  RPK/RP 

*  CALCULATE   DOT  PRODUCT  OF  UNIT  VECTORS 
DOT  =   ((    I*IP   )+(    J*JP   )+(    K*KP    )) 

CALCULATE   FORCES   IN  THE    ' IJK*    SYSTEM 
FHI   =   (MUP/(RP**2))*(R/RP)*(3.  OD+00*DOT*(IP)-(I)) 
FHJ  =   (MTJP/(RP**2))*(R/RP)*(3.  OD+00*DOT*(JP)-(J)) 
FHK  =   (MUP/(RP**2))*(R/RP)*(3.0D+00*DOT*(KP)-(K)) 

*  Transform  FORCES  TO  THE  RSW  SYSTEM 

CALL   IJKRSW(LAN,AL,INC,FHI,FHJ,FHK,FHR,FHS,FHW) 

RETURN 
END 
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SUBROUTINE   SUNPOS(TT,RSI ,RSJ,RSK,RS ,SLAN,SI , SAL, PI) 
THIS   SUBROUTINE  CALCULATES  THE   SUNS   POSITION 

*  VARIABLES  USED  TO  DESCRIBE  THE   SUNS   ORBIT: 
SI      =  SUNS    INCLINATION 

*  SLAN=  SUNS   Longitude  OF  ASCENDING  NODE 
SAP  =   SUNS   ARGUMENT  OF   PERIGEE 

RS      =   SUNS   ORBITAL  RADIUS 
STA  =  SUNS  TRUE  ANOMALY 

*  SAL  =  SUNS  ARGUMENT  OF  LONGITUDE 

DOUBLE   PRECISION     SLAN.SI ,SAL,RS,STA,SAP,TT,RSI ,RSK, 
+       RSJ,RSP,RSQ,RSW,PI 

SI   =  4. 09279709D-01 

SLAN  =  0. OD+00 

SAP  =  0. OD+00 

RS   =    1.4959965D+08 

STA  =  ((2.  0*PI)/(365.0  *  86400.0)   *  TT) 

SAL  =  STA  +  SAP 

*  CALCULATE   SUNS   POSITION   IN    'PQW'    SYSTEM 
RSP  =  RS*DCOS(STA)    • 
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RSQ   =  RS*DSIN(STA) 
RSW  =   0. 0D+00 

*  TRANSFORM   POSITION  TO    'IJK'    SYSTEM 

CALL  PQVIJK(SLAN,SAP,SI,RSP,RSQ,RSW,RSI,RSJ,RSK) 

RETURN 

END 

SUBROUTINE   MONPOS ( TT , RMI , RMJ , RMK , RM , MLAN , MI , MAL , P I ) 
THIS   SUBROUTINE   CALCULATES   THE   MOONS   POSITION 

*  VARIABLES  USED  TO  DESCRIBE  THE   SUNS   ORBIT: 
MI      =  MOONS    INCLINATION 

MLAN=  MOONS   Longitude  OF  ASCENDING  NODE 

MAP  =  MOONS   ARGUMENT  OF  PERIGEE 

RM     =  MOONS   ORBITAL  RADIUS 

MTA  =  MOONS  TRUE  ANOMALY 

MAL  =  MOONS   ARGUMENT  OF  LONGITUDE 

DOUBLE   PRECISION     MI , MLAN, MAL, RM,TM, MTA, RMP, RMQ, RMV, 
+  RMI, RMJ, RMK, MAP, TT, PI 

MI   =  4. 99164166D-01 

RM  =   3. 844D+05 

MLAN  =   0.0 

MTA  =   ((2. 0*PI)/(27. 3   *   3600)    *  TT) 

MAP   =   0. OD+00 

MAL  =  MTA 

CALCULATE   MOON  POSITION   IN    'PQW*    SYSTEM 
RMP   =  RM*DCOS(MTA) 
RMQ  =  RM-DSIN(MTA) 
RMV  =  0 

TRANSFORM  POSITION  TO    'UK*    SYSTEM 

CALL   PQW I JK ( MLAN , MAP , M I , RMP , RMQ , RMV , RM I , RMJ , RMK ) 

RETURN 
END 
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SUBROUTINE   FDRAG( RI , RJ , RK , R , VI , VJ , VK , V , LAN , AP , I , FDR , FDS , FDV , 
+  EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 

+  MM,N,H,HI,HJ,DT) 

THIS   SUBROUTINE   CALCULATES  THE   PERTURBING  FORCE  DUE  TO  DRAG 

THE   FOLLOWING  VARIABLES  ARE   USED  TO  MODEL  THE  ATMOSPHERE: 
RE  =  RADIUS   OF  EARTH 

M  =  MASS   OF   SATELLITE 

=  FRONTAL  SURFACE   AREA  OF   SATELLITE 


AR 

Z 

K 

DENO 
CD 


=   ALTITUDE   OF   SATELLITE 
=   EXPONENTIAL  DECAY  FACTOR 
=   NORMAL  DENSITY 
COEFFICIENT  OF   DRAG 
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DOUBLE  PRECISION  MAG,M,K,FDR,FDS ,FDW,RE , AR , Z ,DENO,CD,DEN, 
+  FDJ,FDK,FDI,RI,RJ,RK,VI,VJ,VK,V,LAN,AP,I,R, 

+  EI,EJ,EK\E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 

+  MM  ,  N  ,  H ,  HI  ,  HJ ,  DT ,  DVR ,  DVS  ,  DVW ,  DVI ,  DVJ ,  DVK 

RE  =  6. 378145D+03 
M  =  1.  OD+02 
AR  =  2.  OD+01 
Z  =  R  -  RE 

DEPENDING  ON  ALTITUDE  SET  ATMOSPHERE  VARIABLES 
IF  (Z.LE.  1.5D+02)  THEN 

K  =  4. 74D-02 

DENO  =  1. 225D+00 

CD  =  1.  OD+00 
ELSEIF  (Z. LE. 5. 5D+02)  THEN 

K  =  3. 4614D-02 

DENO  =  1. 79846D-01 

CD  =  2. OD+00 
ELSE 

K  =  2. 21698D-3 

DENO  =  1. 015484D-07 

CD  =  2. OD+OO 
END  IF 

CALCULATE  ATMOSPHERIC  DENSITY 
DEN  =  DENO  *  DEXP(-K*Z) 

CALCULATE  MAGNITUDE  OF  DRAG  FORCE  AND  LIMIT  IT  TO  1. OE-20 
MAG  =  -(0.  5D+00)*CD*AR*DEN*V*(1.  0D-03)/M 
IF  (ABS(MAG)  .  LT.  1.0D-20)  THEN 

MAG  =  -1.  OD-20 
END  IF 

GIVE  DRAG  FORCE  A  Direction  OF  MINUS  THE  VELOCITY 

FDR  =  0.0 

FDS  =  MAG  *  V 

FDW  =  0.0 

RETURN 

END 

CALCULATE  PERTURBED  NEW  ELEMENTS 


SUBROUTINE  PNEWEL(FR,FS ,FW,H,R, A,E ,N,TA,DT, I , LAN, AL, AP ,P , 

+   MM , MA , E A , TF , T , MU , P I , D I , DA , DE , DMM , DMA , DLAN , DH , DAP ) 
THIS  SUBROUTINE  CALCULATES  THE  NEW  ELEMENTS  FROM  THE  PREVIOUS 
ELEMENTS  ADDED  TO  THE  RATES  OF  CHANGE  FOR  ONE  STEP 

THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

RATE   =  CALCULATES  RATES  OF  CHANGE  OF  ORBITAL  ELEMENTS 

NANGMO  =  NEW  ANGULAR  MOMENTUM  (NEWH) 

NSMA   =  NEW  SEMI -MAJOR  AXIS  (NEWA) 

NECC   =  NEW  ECCENTRICITY  (NEWE) 
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•'; 

NINCL 

=  NEW 

-.'.- 

NASNOD 

=  NEW 

-,V 

NARPER 

=  NEW  . 

-."; 

NMNMO 

=  new  : 

-.'.- 

MEAN'MO 

=  MEAN 

V.- 

NMMAN 

=  NEW  : 

-V 

NEA 

=  new  : 

»v 

NTA 

=  NEW  ' 

v.- 

TFLGHT 

=  TIME 

INCLINATION  (NEWI) 

LONGITUDE  OF  ASCENDING  NODE  (NEWLAN) 

ARGUMENT  OF  PERIGEE  (  NEWAP) 

MEAN  MOTION  (NEWMM) 

MOTION  (MM) 
MEAN  ANOMALY  (NEWMA) 
ECCENTRIC  ANOMALY  (EA) 
TRUE  ANOMALY  (TA) 

OF  FLIGHT  (TF) 

DOUBLE  PRECISION  FR,FS ,FW,DMM,H,R, A,E ,N,TA,DT, I ,LAN,AL, AP ,P, 
+   MM , MA , EA , TF , T , MU , PI , DA , DH , DE , DI , DLAN , DAP , DMA , 
+   NEWH , NEWA , NEWE , NEWI , NEWLAN , NEWAP , NEWMM 

INCREMENT  TIME  OF  FLIGHT  BY  ONE  TIME  STEP  AND  CALCULATE  RATES 
TF  =  TF  +  DT 

CALL  RATES(DH,DA,DE,DI,DLAN,DAP,DMM,DMA,E,MM,R,A,FR,FS,FW, 
+  TA,AL,H,P,T,MU,I) 

CALCULATE  NEW  ELEMENTS 

CALL  NANGMO(H,DT,DH,NEWH) 

CALL  NSMA(A,DT,DA,NEWA) 

CALL  NECC(E,DT,DE,NEWE) 

CALL  NINCL(I,DT,DI,NEWI) 

CALL  NASNODC  LAN , DT , DLAN , NEWLAN ) 

CALL  NARPER(AP,DT, DAP,  NEWAP) 

SET  ELEMENTS  TO  NEW  ELEMENTS 

A  =  NEWA 

l,   —   NEVi  E 

I  =  NEWI 

LAN  =  NEWLAN 

AP  =  NEWAP 

P  =  A  *  (1  -  E ••'•■■-'•- 2) 

MOVE  THE  SATELLITE  ONE  TIME  STEP 

CALL  MEANMO(A,MM,MU) 

CALL  NMNAN ( MA , MM , DT , TF , DMA ,  P I ) 

CALL  NEA(MA,E,EA) 

CALL  NTA(EA,E,TA,PI) 

CALL  TFLGHT(MM,MA,TF) 

AL  =  TA  +  AP 

RETURN 

END 
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CALCULATE  THE  RATES  OF  CHANGE  OF  THE  ORBITAL  Elements 


SUBROUTINE  RATES( DH,DA,DE ,DI , DLAN, DAP, DMM, DMA, E ,MM,R, A,FR,FS ,FW, 
+  TA,AL,H,P,T,MU,I) 
THIS  SUBROUTINE  Calls  THE  FOLLOWING  SUBROUTINES  TO  CALCULATE  THE 
TIME  RATE -OF-  CHANGE  OF  THE  ORBITAL  ELEMENTS: 
RSMAX  =  RATE-OF-CHANGE  OF  THE  SEMI -MAJOR  AXIS  (DA) 
RECC   =  RATE -OF -CHANGE  OF  THE  ECCENTRICITY  (DE) 
RLNC    =  RATE-OF-CHANGE  OF  THE  INCLINATION  (DI) 
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RLAN   =  RATE-OF-CHANGE  OF  THE  Longitude  OF  THE  ASCENDING  NODE 

(DLAN) 
RAP    =  RATE-OF-CHANGE  OF  THE  ARGUMENT  OF  PERIGEE  (DAP) 
RMM    =  RATE -OF -Change  OF  THE  MEAN  MOTION  (DMM) 
RMA    =  RATE-OF-CHANGE  OF  THE  MEAN  ANOMALY  (DMA) 
RANGMO  =  RATE-OF-CHANGE  OF  THE  ANGULAR  MOMENTUM  (DH) 

DOUBLE  PRECISION  DH,DA,DE ,DI , DLAN, DAP , DMM, DMA, E ,MM,R, A,FR,FS ,FW, 
+  TA,AL,H,P,T,MU,I 

CALL  RSMAX( E , MM , R , A , FR , FS , DA , TA) 

CALL  RECC(E,MM,R,A,FR,FS,TA,DE) 

CALL  RINC(E,MM,R,A,FW,AL,DI) 

CALL  RLAN( E , MM , R , A , I , FW , AL , DLAN) 

CALL  RAP(E,MM,R,A,I,H,P,AL,TA,FR,FS,FW,DAP) 

CALL  RMM(MM,A,DMM,DA,MU) 

CALL  RMA(E,MM,R,A,TA,DMM,FR,FS,DMA,T) 

CALL  RANGMO (R,FS,FW,DH) 

RETURN 

END 

SUBROUTINE  RANGMO( R , FS , FW , DH ) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE 

ANGULAR  MOMENTUM 

DOUBLE  PRECISION  FS ,FW,DHW,DHS ,DH ,R 

DHW  =  R  *  FS 

DHS  =  R  *  FW 

DH  =  DSQRT((DHW**2)  +  (DHS **2)) 

RETURN 

END 

SUBROUTINE  RSMAX(E ,MM ,R , A ,FR,FS ,DA,TA) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  SEMI -MAJOR 

AXIS 

DOUBLE  PRECISION  DA,FR,FS ,E,MM,R,A,TA,ET 

*     TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (E.  GT.  0.  9)  THEN 

ET  =  0.  9 
ELSE 

ET  =  E 
END  IF 
DA  =  ((2.0D+00*E  -DSIN(TA))/(MM  *DSQRT( 1.  0D+00-(ET**2) ) ) )*FR  + 
+    ( (  2.  0D+00*A*DSQRT(  1.  0D+00-(E  **2)  )  )/(MM  *R)  )*FS 
RETURN 
END 
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SUBROUTINE   RECC ( E , MM , R , A , FR , FS , TA , DE ) 

THIS   SUBROUTINE   CALCULATES  THE  RATE   OF  CHANGE   OF  THE  ECCENTRICITY 

DOUBLE   PRECISION      DE ,FR ,FS ,E ,MM,R, A ,TA,ET 

TRAP   (E)    SO  DENOMINATOR  DOES  NOT  GOTO   ZERO 
IF    (E.  LT.  0.  1)    THEN 

ET  =   0.  1 
ELSE 

ET  =  E 
END  IF 
DE  =   ((DSQRTQ.OD+00    -    (E  **2) )*SIN(TA) )/(MM  *A))*FR  + 
+        ((DSQRT(1.0D+00    -    (E  **2) ) )/(MM  *ET*( A** 2) ) )* 
+        ((A**2)*(1.0D+00    -    (E   **2))/(R)    -    (R))*FS 
RETURN 
END 

SUBROUTINE  RLAN( E , MM ,R, A , I ,FW, AL, DLAN) 

THIS   SUBROUTINE   CALCULATES  THE  RATE   OF  CHANGE   OF  THE   LONGITUDE 

OF  THE  ASCENDING  NODE 

DOUBLE   PRECISION     DLAN,FW,E ,MM,R, A, I , AL,ET, IT 

TRAP   (El    AND   (I)    SO  DENOMINATOR  DOES   NOT  GOTO   ZERO 
IF    (E.GT.  0.  9)    THEN 

ET  =   0.  9 
ELSE 

ET  =  E 
END  IF 
IF   (I.  LT.  0.  01 745)    THEN 

IT  =   0. 01745 
ELSE 

IT  =    I 
ENLIF 

DLAN   =   (R*FW*DSIN(AL))/(MM  *( A**2)*DSQRT( 1.  OD+00    -    (ET**2))* 
+        DSIN(IT)) 
RETURN- 
END 
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SUBROUTINE  RAP(E ,MM,R, A, I ,H,P, AL,TA,FR,FS ,FW,DAP) 

THIS    SUBROUTINE   CALCULATES   THE   RATE   OF   CHANGE   OF  THE   ARGUMENT 

OF   PERIGEE 

DOUBLE   PRECISION     DAPR.DAPS ,DAPW,DAP ,FR,FS ,FW,E ,MM,R, I ,H,P, AL,TA, 

+       ET,A,IT 

TRAP   (E)    AND   (I)    SO  DENOMINATOR  DOES   NOT  GOTO   ZERO 
IF   (I.LT.  0.  01745)   THEN 
IT  =   0. 01745 

ELSE 

IT  =    I 

en::.' 


ORB17300 
ORB17310 
ORB17320 
ORB17330 
ORB17340 
ORB17350 
ORB17360 
ORB17370 
ORB17380 
ORB17390 
ORB17400 
0RB17410 
ORB17420 
ORB17430 
ORB17440 
ORB17450 
ORB17460 
0RB17470 
ORB17480 
ORB17490 
ORB17500 
ORB17510 
ORB17520 
ORB17530 
ORB17540 
ORB17550 
ORB17560 
ORB17570 
ORB17580 
ORB17590 
ORB17600 
ORB17610 
ORB17620 
ORB17630 
ORB17640 
ORB17650 
ORB17660 
ORB17670 
ORB17680 
ORB17690 
ORB17700 
ORB17710 
ORB17720 
ORB17730 
ORB17740 
ORB17750 
ORB17760 
ORB17770 
ORB17780 
ORB17790 
ORB17800 
ORB17810 
ORB17820 
ORB17830 
ORB  17840 
ORB17850 


57 


IF  (E.  GT.  0.  9)  THEN 

ET  =  0.  9 
ELSEIF  (E.  LT.  0.  1)  THEN 

ET  =  0.  1 
ELSE 

ET  =  E 
END  IF 

DAPR  =  (-DSQRT(  1.0+00  -  (E  **2) )*DCOS(TA) )/(MM  *A*ET)  *  FR 
DAPS  =  (P/(ET*H))*(DSIN(TA))* 
+    (1.0D+00  +  1.  0D+00/(l.  0D+00  +  ET*DCOS(TA) ) )  *FS 

DAPW  =  (-R*(l. OD+00/DTAN(IT))*DSIN(AL))/ 
+      (MM  *(A**2)*DSQRT(1.  0D+00  -  (ET**2)  )  )*FW 
DAP  =  DAPR  +  DAPS  +  DAPW 
RETURN 
END 

SUBROUTINE  RINC(E,MM,R, A,FW, AL.DI) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  INCLINATION 

DOUBLE  PRECISION  DI ,FW,E,MM,R,A,AL,ET 

*     TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 
IF  (E.GT.  0.  9)  THEN 

ET  =  0.  9 
ELSE 

ET  =  E 
END  IF 

DI  =  (R*FW*DCOS(AL))/(MM  *(A**2)*DSQRT( 1. OD+00  -  (ET**2))) 
RETURN 
END 

SUBROUTINE  RMM( MM , A , DMM , DA , MU) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  MEAN  MOTION 

DOUBLE  PRECISION  DMM, DA, MM, A, MU 

DMM  =r(-3. OD+00*MU)/(2. 0D+00*MM  *(A**4)))*  DA 

RETURN 
END 

SUBROUTINE  RMA( E , MM , R , A , TA , DMM , FR , FS , DMA , T) 

THIS  SUBROUTINE  CALCULATES  THE  RATE  OF  CHANGE  OF  THE  MEAN  Anomaly 

DOUBLE  PRECISION  DMAA.DMAB ,DMAC ,DMAD ,DMM,FR,FS ,DMA,E ,MM,R, A,TA, 

+   ET,T 

TRAP  (E)  SO  DENOMINATOR  DOES  NOT  GOTO  ZERO 

IF  (E.  GT.  0.  9)  THEN 

ET  =  0.  9 
ELSEIF  (E.  LT.  0.  1)  THEN 
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ET  =  0.  1 

ELSE 

ET  =  E 

END  IF 

DMA  =  (-1.  0D+00/(MM  *A))* 
+       (((2. 0D+00*R)/A)  -  ((1  -  (E  **2))/ET)*DCOS(TA))  *  FR  - 
+     (1-(E  **2))/(MM  *A*ET)*(1+  R/(A*( 1-(E**2))))*(SIN(TA)*FS) 
+       (T  *  DMM) 

RETURN 

END 
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*  CALCULATE  THE  NEW  ORBITAL  ELEMENTS 

SUBROUTINE  NSMA( A,DT,DA,NEWA) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  SEMI -MAJOR  AXIS 

DOUBLE  PRECISION  DA , DT , A , NEVA 

NEVA  =  A  +  DA*DT 

RETURN 

END 

Ay?Ay«vvr?vyf  Av«v«vyvVTyryryry*y*y?Ayryyyr^ 

SUBROUTINE  NECC( E , DT , DE , NEVE ) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  ECCENTRICITY 

DOUBLE  PRECISION  DE,DT,E,NEVE 

NEVE  =  E  +  DE-DT 

RETURN 

END 

SUBROUTINE  NINCL( I ,DT,DI ,NEVI ) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  INCLINATION 

DOUBLE  PRECISION  DI,DT,I,NEVI 

NEW I  =  I  +  DI*DT 
RETURN- 
END 

iVyryfyryriViVVwryryry.-VwVywWrVryry.-yryryr^ 

SUBROUTINE  NASNOD( LAN , DT , DLAN , NEVLAN) 

*  THIS  SUBROUTINE  CALCULATES  THE  NEW  LONGITUDE  OF  THE  ASCENDING 
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D0U3LE  PRECISION  DLAN, DT, LAN, NEWLAN 

NEVLAN  =  LAN  +  DLAN-DT 

RETURN 

END 
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SUBROUTINE  NARPER( AP , DT , DAP , NEWAP ) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  ARGUMENT  OF  PERIGEE 

DOUBLE  PRECISION  DAP ,DT,AP , NEWAP 

NEWAP  =  AP  +  DAP*DT 

RETURN 

END 

*iV*Vfyryfyr7VyryrycyciV-V*yrycyr*ycyWryrVr^^^ 

SUBROUTINE  NMNAN(MA,MM,DT,TF,DMA,PI) 
*     THIS  SUBROUTINE  CALCULATES  THE  NEW  MEAN  Anomaly 

DOUBLE  PRECISION  DMM ,FR,FS ,DMA,DT,MA,E ,R, A ,TA,MM,TF,T, PI 

MA  =  MM*(TF)   +  DMA*DT 
IF  (MA  . GT.  (2*PI))  THEN 
MA  =  MA  -  (2*PI) 

ENDIF 

RETURN 
END 

SUBROUTINE  NMNM0( MM , DMM , DT , NEWMM ) 

THIS  SUBROUTINE  CALCULATE  THE  NEW  MEAN  MOTION 

DOUBLE  PRECISION  DMM , DT , MM , NEWMM 

NEWMM  =  MM  +  DMM*DT 

RETURN 

END 
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SUBROUTINE  NEA(MA,E,EA) 

THIS  SUBROUTINE  CALCULATES  THE  NEW  ECCENTRIC  ANOMOLY  BY  USING 

NEWTONS  METHOD  OF  ROOT  FINDING 

DOUBLE  PRECISION  EAN,MAN,MA,E ,EA,DIFF 

•     LET  (EA)  EQUAL  (MA)  FOR  INITIAL  GUESS  AT  ROOT 
EA  =  MA 

EAN  =  EA  +  (MA  -  EA  +E*DSIN(EA) )/( 1.  OD+00  -  E*DCOS(EA)) 
MAN  =  EAN  -  E-SIN(EAN) 

CHECK  DIFFERENCE  (DIFF) 
DIFF  =  ABS(MA  -MAN) 
EA  =  EAN 

CONTINUE  TO  INTERATE  UNTIL  DIFFERENCE  IS  NEGLIGIBLE 
200   IFCDIFF. GT. 0. 00000C0001)  THEN 
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EAN  =  EA  +  (MA  -  EA  +  E*DSIN(EA) )/( 1.  OD+00  -  E*DCOS(EA))  ORB19540 

MAN  =  EAN  -  E--'--DSIN(EAN)  ORB  19550 

EA  =  EAN  0RB19560 

DIFF  =  ABSCMA  -  MAN)  ORB19570 

GOTO  200  ORB  19580 

ENDIF  ORB19590 

EA  =  EAN  ORB  19600 

RETURN  ORB  196 10 

END  ORB19620 

ORB19630 

SMnV**iWrV«v****}V********tov**toV^^  ORB19640 

ORB19650 
SUBROUTINE  NTA(EA,E ,TA,PI)  ORB  19660 

THIS  SUBROUTINE  CALCULATES  THE  NEW  TRUE  Anomaly  ORB19670 

ORB19680 
DOUBLE  PRECISION  EA,E,TA,PI  ORB19690 

ORB19700 
TA  =  DACOS((E  -  DCOS(EA))/(E*DCOS(EA)  -  1.  OD+00))  OR319710 

IF  (EA.  GT.  PI)  THEN  ORB19720 

TA  =  (2--'PI)  -  TA  ORB  19730 

ENDIF  ORB19740 

RETURN  ORB 19 750 

END  ORB19760 

ORB19770 
****»v**yc*>vvc********'fofr**fov*****v^  ORB  19  780 

ORB19790 

SUBROUTINE  NANGMO(H ,DT,DH,NEWH)  ORB19800 

THIS  SUBROUTINE  CALCULATES  THE  NEW  ANGULAR  MOMENTUM  ORB  198 10 

ORB19820 
DOUBLE  PRECISION  DH,DT,H,NEWH  ORB  19830 

0RB19840 
NEWH  =  H  +  DH-DT  ORB  19 850 

RETURN  ORB  19860 

END  ORB  19870 

ORB19880 

ORB19900 

SUBROUTINE  INTSUM(TFEA ,TFSU,TFMO,TFDRA,TDI ,TDA,TDE ,TDMM,TDMA,  ORB19910 
+  TDLAN,TDH,TDAP)  ORB  19920 

THIS  SUBROUTINE  INITIALIZES  THE  SUMS  OF  FORCES  AND  ELEMENT  CHANGESORB 19930 

ORB19940 
TFEA  =  0.  0  ORB19950 

TFSU  =  0.  0  ORB19960 

TFMO  =  0.  0  ORB19970 

TFDRA  =  0.  0  ORB19980 

TDI  =  0.  0  ORB19990 

TDA  =  0.  0  ORB20000 

TDE  =  0. 0  ORB20010 

TDMM  =  0.0  ORB20020 

TDMA  =0.0  ORB20030 

TDLAN  =0.0  ORB20040 

TDH  =  0.0  ORB20050 

TDAP  =0.0  ORB20060 

RETURN  ORB20070 

END  ORB20080 


-)--.*-»;-*• -*•-.'- 
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■••■■  CALCULATE  THE  NEW  POSITION  AND  VELOCITY  VECTORS 

SUBROUTINE  NPOS(RI ,RJ,RK,R,LAN,AP , INC ,  TA,A,E) 
THIS  SUBROUTINE  CALCULATES  THE  NEW  POSITION  VECTOR 

DOUBLE  PRECISION  XW, YW,ZW, INC.RI ,RJ,RK,R,LAN, AP,TA,A,E 

CALCULATE  POSITION  VECTOR  IN  'PQW'  SYSTEM 
R  =  (A*(l  -  (E**2)))/(l  +  E*DCOS(TA)) 
XW  =  R--'-DCOS(TA) 
YW  =  R*DSIN(TA) 

ZW  =  0 

TRANSFORM  POSITION  TO  'UK'  SYSTEM 

CALL  PQWIJK(LAN,AP,INC,XW,YW,ZW,RI,RJ,RK) 

R  =  DSQRT((RI**2)  +  (RJ**2)  +  (RK**2)) 

RETURN 
END 

SUBROUTINE  NVEL(E ,P,TA,LAN,AP, INC, VI , VJ,VK,V,MU) 
THIS  SUBROUTINE  CALCULATES  THE  NEW  VELOCITY  VECTOR 

DOUBLE  PRECISION   INC , VP, VQ, VW,MU,E , P,TA, LAN, AP , VI , VJ, VK, V 

CALCULATE  VELOCITY  IN  'PQW'  SYSTEM 
VP  =  DSQRT(MU/P)*(-DSIN(TA)) 
VQ  =  DSQRT(MU/P)*(E  +  DCOS(TA)) 
VW  =  0.  OD+00 

*  TRANSFORM  VELOCITY  INTO  '  UK'  SYSTEM 
CALL  PQWI JK( LAN , AP , INC , VP , VQ , VW , VI , VJ , VK) 
V  =  DSQRT((VI**2)  +  (VJ**2)  +(VK**2)) 

RETURN 
END 

VELOCITY  CHANGE 

SUBROUTINE  CHGVEL(DT,PER, AL,LAN, AP, I ,RI ,RJ,RK,R, 
+   VI,VJ,VK,V,MU,PI,H,A,E,N,TA,P,MM,MA,EA, 

+   TF,T,NUM,RIRAY,RJRAY,RKRAY,RARAY,TARAY,AINRAY,APRAY,TIMRAY, 
+   TT,EI,EJ,EK,LP,HI,HJ,I0PT1,TFEA,TFSU,TFM0,TFDRA, 
+    TD I , TDA , TDE , TDMM , TDM A , TDLAN , TDH , TDAP ) 

*  THIS  SUBROUTINE  CALCULATE  VELOCITY  CHANGES 

*  THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

TACHG  =  RETURNS  TRUE  ANOMALY  FOR  VELOCITY  CHANGE  LOCATION  (CHTA) 

AND  AN  INDICATOR  OF  LOCATION  ( ITA) 
CALCEL  =  CALCULATE  Orbital  ELEMENTS 
UNPRET  =  CALCULATE  UNPERTURBED  ORBIT 


ORB20090 
0RB20100 
ORB20110 
0RB20120 
ORB20130 
0RB20140 
ORB20150 
ORB20160 
ORB20170 
0RB20180 
ORB20190 
ORB20200 
ORB20210 
ORB20220 
ORB20230 
ORB20240 
ORB20250 
ORB20260 
ORB20270 
ORB20280 
ORB20290 
ORB20300 
ORB20310 
ORB20320 
ORB20330 
ORB20340 
ORB20350 
ORB20360 
ORB20370 
ORB20380 
0RB20390 
ORB20400 
ORB20410 
ORB20420 
ORB20430 
ORB20440 
ORB20450 
ORB20460 
ORB20470 
0RB20480 
ORB20490 
ORB20500 
ORB20510 
ORB20520 
0RB20530 
ORB20540 
ORB20550 
ORB20560 
0RB20570 
ORB20580 
ORB20590 
ORB20600 
ORB20610 
ORB20620 
ORB20630 
ORB20640 


NPOS  =  CALCULATE  NEW  POSITION 

NVEL  =  CALCULATE  NEW  VELOCITY 

•  STORE  =  STORE  POSITION  AND  ELEMENTS  IN  ARRAYS 
ENERGY  =  ENERGY  OF  SATELLITE 

ECC  =  ECCENTRICITY 

SMAXIS  =  SEMI -MAJOR  AXIS 

DOUBLE  PRECISION  T,DT,PER, AL,LAN, AP, I ,RI ,RJ,RK ,R , VI , VJ, VK, V, 
+  MU,PI,H,A,E,N,TA,P,MM,MA,EA,TF,TT, 

+   NEWVI,NEWVJ,NEWVK,NEWV,VMAX,CHTA,EI,EJ,EK,LP,HI,HJ,VCIR, 
+   DI , DE , DA , DMM , DMA , DLAN , DH , DAP , NEWE I , NEWE J , NEWEK , NEWE , NEWENR , 
+   NEWA,NEWRP,R£ 

DIMENSION  RARAY(500) ,TARAY(500) ,RIRAY(500) ,RJRAY(500) , 
+  RKRAY(500) ,AINRAY(500) ,APRAY(500) ,TIMRAY(500) 

CHARACTER* 1 , YORN , PYORN 

RE  =  6.  3782D+03 

PROMPT  THE  USER  FOR  THE  VELOCITY  Change  LOCATION 
CALL  TACNG(PI,CHTA,ITA) 

SET  TIME  COUNTER  TO  ONE  TIME  STEP 
T  =  DT 

■     ROTATE  TO  THE  VELOCITY  CHANGE  LOCATION 

THIS  IS  IDENTICAL  TO  THE  Unperturbed  ORBIT  WITH  THE  EXCEPTION 
THAT  A  COMPLETE  ORBIT  IS  NOT  CALCULATED 
PRINT*, 'ROTATE  TO  VELOCITY  CHANGE  LOCATION ' 
IF  ((ITA.  EQ.  2)  .OR.  ( ITA.  EQ.  3))  THEN 
PRINT* /BEFORE  TA  ='  ,TA 
IF  (TA  .  GT.  6.  21)  THEN 

TA  =  TA  -  (2*PI) 
END  IF 
250     IF((T.  LE.  PER).  AND.  (TA.  LT.  CHTA))  THEN 

•  PRINT*,' TA  =' ,TA 
NUM  =  NUM  +  1 

TT  =  TT  +  DT 

CALL  NEWELT(MM,MA,E,EA,TA,TF,DT,PI,PER) 
CALL  NPOS(RI,RJ,RK,R,LAN,AP,I,  TA,A,E) 
CALL  NVEL(E,P,TA,LAN,AP,I,VI,VJ,VK,V,MU) 
CALL  STORE (RI,RJ,RK,R,TA,RIRAY, R JRA Y , RKRAY , RARA Y , 
+  TARAY,NUM,I,AP,AINRAY,APRAY,TT,TIMRAY) 

T  =  T  +  DT 
GOTO  250 
END  IF 
IF  (TF  . GE.  PER)  THEN 

TF  =  TF  -  PER 
ENDIF 
END  IF 

PRINT  ESCAPE  VELOCITY  AND  CIRCULAR  VELOCITY  FOR  Reference 

CALL  EXCMS( 'CLRSCRN' ) 

PRINT*, 'AFTER  TA  =' ,TA 

PRINT*, 'THIS  SHOULD  BE  THE  DESIRED  RADIUS  RP  OR  RA' 
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160   PRINT--"-, 'RADIUS  ='  ,R 

PRINT--'-, 'VELOCITY  ='  ,V 

VMAX  =  DSQRT(2. 0*(MU  /  R)) 

PRINT-", 'MAX  VELOCITY  AT  THIS  RADIUS  IS:  '  ,VMAX 

VCIR  =  DSQRTCMU/R) 

PRINT"--, 'CIRCULAR  VELOCITY  AT  THIS  RADIUS  IS  :  '  ,VCIR 

PROMPT  USER  TO  CHANGE  VELOCITY  IN  ORBITAL  PLANE 

PRINT*, 'DO  YOU  WANT  TO  CHANGE  THE  VELOCITY  IN  THE  ORBITAL  PLANE?' 

PRINT*, 'ENTER  "Y"  OR  "N"  :  ' 

READ* , PYORN 

PRINT'S  PYORN 

IF  (PYORN  .  EQ.  'Y')  THEN 

PRINT*, 'GIVE  THE  TOTAL  CHANGE  IN  VELOCITY,  I.E.  5.0  KM.  ' 

PRINT*, 'THE  PROGRAM  WILL  FIGURE  OUT  THE  FINAL  VELOCITY  VECTOR' 

PRINTS, '   ENTER  VELOCITY  CHANGE: ' 

READ*,CHGV 

PRINT*,  CHGV 


IN  THE  ORBITAL  PLANE 


CALCULATE  NEW  VELOCITY  FOR  CHANGE 
NEWVI  =  VI  +  (CHGV  *  VI  /  V) 
NEWVJ  =  VJ  +  (CHGV  *  VJ  /  V) 
NEWVK  =  VK  +  (CHGV  *  VK  /  V) 


Velocity  CHANGE  OUT  OF  ORBITAL  PLANE 
ELSE IF  (PYORN  .  EQ.  'N')  THEN 

PRINT*,'  ENTER  THE  NEW  VELOCITY  VECTOR: ' 

PRINT*,'  ENTER  THE  NEW  VI' 

RE AD*, NEW VI 

PR  I  NT*,  NEW  VI 

PRINT*,'  ENTER  THE  NEW  VJ' 

READ* , NEWVJ 

PRINT'-,  NEWVJ 

PRINT'S'  ENTER  THE  NEW  VK* 

RE AD*, NEWVK 

PRINT*, NEWVK 

NUM  =  1 

ITA  =  3 
ELSE 

CALL  EXCMS( 'CLRSCRN' ) 

GOTO  260 
END  IF 

PRINT  NEW  VELOCITY  FOR  USER  TO  CHECK 

NEWV  =  DSQRT((NEWVI**2)  +  (NEWVJ**2)  +  (NEWVK**2)) 

PR I NT* , ' NEW  V I  = ' , NEWV I 

PRINT'S 'NEW  VJ  =',NEWVJ 

PRINT*, 'NEW  VK  =' , NEWVK 

PRINT'S 'NEW  V  ='  ,NEWV 

PRINT*,'  ARE  THESE  VALUES  THE  ONES  YOU  WANT?' 

PRINT*, 'ENTER  "Y"  OR  "N"  :' 

READ* , YORN 

PRINT*, YORN 

IF  (YORN   .EQ.  'N')  THEN 

CALL  EXCMSC  CLRSCRN') 

GOTO  260 
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* 


END  IF 

CHECK  FOR  VALID  VELOCITY 
IF  (  NEWV  . GT.  VMAX)  THEN 

PRINTS, 'YOUR  VELOCITY  IS  TO  GREAT  !! ' 

GOTO  260 
END  IF 

Calculate  PERIGEE  RADIUS  TO  SEE  IF  SATELLITE  WILL  IMPACT  EARTH 

CALL  ENERGY( NEWV , R , MU , NEWENR ) 

CALL  ECC ( RI , R J , RK , R , NEWVI , NEWV J , NEWVK , NEWV , NEWE I , NEWE J , NEWER , 
+         NEWE,MU) 

CALL  SMAXIS(MU, NEWENR, NEWA) 

NEWRP  =  NEW  A*  (1.0  -  NEWE) 

IF  (NEWRP  .  LE.  RE)  THEN 

PRINT*, 'YOUR  VELOCITY  AT  THIS  POINT  IS  TO  SMALL!!! ' 

PRINT*, 'THE  SATELLITE  WILL  IMPACT  THE  EARTH!!  ' 

PRINT*,' THE  SATELLITES  RADIUS  OF  PERIGEE  WOULD  BE  ', NEWRP 

PRINT*, 'A  NEW  VELOCITY  WILL  HAVE  TO  BE  ENTERED!! ' 

GOTO  260 

END  IF 

ACCEPT  NEW  VELOCITY 

VI  =  NEWVI 
VJ  =  NEWVJ 
VK  =  NEWVK 
V  =  NEWV 

CALCULATE  NEW  ELEMENT  WITH  NEW  VELOCITY  AND  SET  TIME  STEP 
CALL  CALCEL(  RI ,RJ ,RK ,R, VI , VJ , VK, V ,EI ,EJ ,EK,E , A , I ,LAN, L? ,TA , 
+  PER,EA,MA,AP,AL,TF,P,PI,MU,MM,N,H,HI,HJ) 

DT  =  PER/50.  0 
T  =  DT 

THE  FOUR  Different  CASES  OF  VELOCITY  CHANGES  FOLLOWS: 

VELOCITY  CHANGE  AT  PERIGEE,  AND  NEWV  >  V  Circular 

IF((ITA.  EQ.  1).  AND.  (NEWV.  GT.  VCIR)  )THEN 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA,EA,TF,T,NUM,RIRAY,RJRAY,RKRAY,RARAY, 

+  TARAY,AINRAY,APRAY,TIMRAY,TT) 

Change  VELOCITY  AT  PERIGEE,  AND  NEWV  <=  V  CIRCULAR 

APOGEE  AND  PERIGEE  SWAP 

ELSEIF  ((ITA.  EQ.  1).  AND.  (NEWV.  LE.  VCIR)  )THEN 

CLEAR  PREVIOUS  PLOTS 
NUM  =  1 

CALL  STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY,TARAY, 
+  NUM ,  I ,  AP ,  AINRAY ,  APRAY , TT , TIMRAY) 

T  =  PER/2  , 

STEP  SATELLITE  TO  NEW  PERIGEE,  ONLY  A  HALF  ORBIT 
CALL  UNPRET( DT , PER , AL , LAN , AP , I , RI , R J , RK , R , VI , V J , VK , V , 
+  MU,PI,H,A,E,N,TA,P,MM, 
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MA , E A , TF , T , NUM , R I RA Y , R JRAY , RKRAY , RARAY , 
TARAY .AINRAY, APRAY, TIMRAY, TT) 


RESET  TIME  COUNTER  TO  ONE  TIME  STEP 
T  =  DT 

CALCULATE  COMPLETE  NEXT  ORBIT 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA , E A , TF , T , NUM , R IRAY , R JRAY , RKRAY , RARAY , 

+  TARAY,AINRAY,APRAY,TIMRAY,TT) 

CHANGE  VELOCITY  AT  APOGEE,  AND  NEW  V  <  V  CIRCULAR 
ELSEIF  ((ITA.EQ.  2)  .AND.  (NEWV  .  LT.  VCIR))  THEN 

T  =  PER/ 2 

FINISH  ORBIT 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA , E A , TF , T , NUM , R I RAY , R JRAY , RKRAY , RARAY , 

+  TARAY,AINRAY,APRAY,TIMRAY,TT) 

CHANGE  VELOCITY  AT  Apogee,  AND  NEWV  >=  V  CIRCULAR 

OR  AT  ANY  OTHER  TRUE  Anomaly 

ELSEIF  (((ITA.EQ.  2).  AND.  (NEWV.  GE.  VCIR))  .OR.  (ITA.EQ.  3))  THEN 

IF  (TA  .  GT.  6.  21)  THEN 
TA  =  TA  -  (2*PI) 

END  IF 

CLEAR  PREVIOUS  ORBITS  AND  STEP  SATELLITE  TO  NEW  PERIGEE 

T  =  TF 

NUM  =  1 

CALL  STORE ( R I , R J , RK , R , TA , RIRAY , R JRAY , RKRAY , RARAY , TARAY , 
+  NUM, I, AP,AINRAY, APRAY, TT, TIMRAY) 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA , EA , TF , T , NUM , RIRAY , R JRAY , RKRAY , RARAY , 

+  TARAY,AINRAY,APRAY,TIMRAY,TT) 

IF  (TF  .GE.  PER)  THEN 
TF  =  TF  -  PER 

END  IF 

CALCULATE  COMPLETE  NEXT  ORBIT 
T  =  DT 

CALL  UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VIJVJ,VK,V, 
+  MU,PI,H,A,E,N,TA,P,MM, 

+  MA , E A , TF , T , NUM , R I RA Y , R JRA Y , RKRAY , RARAY , 

+  TARAY,AINRAY,APRAY,TIMRAY,TT) 

END  IF 
RETURN 
END 
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SUBROUTINE  TACNG( PI , CHTA , ITA) 

THIS  SUBROUTINE  Asks  THE  USER  FOR  VELOCITY  CHANGE  LOCATION 

DOUBLE  FRECISION  CHTA, PI 

CALL  EXCMS('CLRSCRN') 


WHERE  DO  YOU  WANT  TO  CHANGE  THE  VELOCITY?' 

1.  AT  CURRENT  PERIGEE' 

2.  AT  CURRENT  Apogee' 

3.  AT  A  SPECIFIC  TRUE  Anomaly' 


ENTER  "1",  "2"  OR  "3"' 


PRINT*, 
PRINT*, 

PRINT*, 
PRINT*, 
PRINT*, 
READ*, ITA 
PRINT*, ITA 

*  SET  TRUE  ANOMALY  CHANGE  LOCATION  (CHTA)  TO  DESIRED  LOCATION 
IF  (ITA  .EQ.  1)  THEN 

CHTA  =  0.0 
END  IF 
IF  (ITA  .EQ.  2)  THEN 

CHTA  =  PI 
END  IF 
IF  (ITA  .EQ.  3)  THEN 

PRINT*, 'AT  WHAT  TRUE  ANOMALY  DO  YOU  WANT  TO  CHANGE  THE' 

PRINT*, 'VELOCITY?' 

PRINT*,' ENTER  TRUE  ANOMALY  IN  DEGREES' 

RE AD*, CHTA 

PR I NT*, CHTA 

CHTA  =  CHTA  *  PI  /  180 
END  IF 
RETURN 
END 

*  OUTPUT  PLOTS 
yryfyr?v*y:y.-^\-y,-yry.-yryrVw?yfyry^ 

SUBROUTINE  PLOTS (RIRAY,RJRAY,RKRAY,RARAY,TARAY,NUM, PI , INC, LP, A, 
+  E,TF,AINRAY,APRAY,TIMRAY,TFEA,TFSU,TFMO,TFDRA, 

+  PER , TD I , TD A , TDE , TDMM , TDMA , TDLAN , TDH , TD AP , 

+  MM,MA,LAN,H,AP,R,V) 

THIS  SUBROUTINE  ASKS  THE  USER  FOR  THE  TYPE  OF  OUTPUT  THAT  IS 

*  DESIRED  PERIFOCAL,  GROUND  TRACK  OR  TO  SKIP  THE  PLOT. 

THE  FOLLOWING  SUBROUTINES  ARE  CALLED: 

*  PERIF  =  PLOT  PERIFOCAL  ORBIT 
GRTRK  =  PLOT  GROUND  TRACK 
DATE   =  DISPLAYS  DATA  ON  PLOT 
TEC618  =  SET  Disspla  TO  TEC  618  OUTPUT 

*  ENDPL  =  END  THIS  DISSPLA  PLOT 

REFER  TO  DISSPLA  USER"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 

*  SUBROUTINES 

DOUBLE  PRECISION  PI , A , E , INC ,LP ,TF ,FER , MM, MA, LAN, H , AP ,R, V 
DIMENSION  RIRAY(500),RJRAY(500),RKRAY(500),RARAY(500),TARAY(500), 


ORB22890 
0RB22900 
0RB22910 
ORB22920 
ORB22930 
ORB22940 
ORB22950 
ORB22960 
ORB22970 
ORB22980 
ORB22990 
ORB23000 
ORB23010 
ORB23020 
ORB23030 
ORB23040 
ORB23050 
0RB23060 
ORB23070 
ORB23080 
0RB23090 
0RB23100 
ORB23110 
ORB23120 
ORB23130 
ORB23140 
ORB23150 
ORB23160 
ORB23170 
ORB23180 
ORB23190 
ORB23200 
ORB23210 
ORB23220 
0RB23230 
ORB23240 
ORB23250 
ORB23260 
ORB23270 
ORB23280 
ORB23290 
0RB23300 
ORB23310 
ORB23320 
ORB23330 
ORB23340 
0RB23350 
0RB23360 
0RB23370 
ORB23380 
ORB23390 
ORB23400 
0RB23410 
0RB23420 
ORB23430 
ORB23440 
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+     AINRAY(500) ,APRAY(500) ,TIMRAY(500) 
CHARACTER* 1,Y0RN 

CALL  EXCMS( 'CLRSCRN') 

CALCULATE  SINGLE  PRECISION  VARIABLES 
SPI  =  SNGL(PI) 
SA  =  SNGL(A) 
SE  =  SNGL(E) 
SINC  =  SNGL(INC) 
SLP  =  SNGL(LP) 
STF  =  SNGL(TF) 
SPER  =  SNGL(PER) 
SMM  =  SNGL(MM) 
SMA  =  SNGL(MA) 
SLAN  =  SNGL(LAN) 
SH  =  SNGL(H) 
SAP  =  SNGL(AP) 
SV  =  SNGL(V) 
SR  =  SNGL(R) 


340 


350 


* 


PROMPT  USER  FOR  DISPLAY  TYPE 


PRINT*, 

PRINT*, 
PRINT*, 
PRINT*, 
PRINT*, 


WHAT  TYPE  OF  Display  IS  DESIRED: ' 

1.  PERIFOCAL* 

2.  GROUND  TRACK' 

3.  SKIP  PLOT' 
ENTER  1,2,3,4:  ' 


READ*, INPUT 
PRINT350, INPUT 
F0RMAT(I4) 

CALL  TEK618 

CALL  APPROPRIATE  PLOT 
IF  (INPUT  . EQ.  1)  THEN 

CALL  PERIF( RARAY , TARAY , NUM , SPI , SINC , SLP , SA , SE ) 
ELSEIF  (INPUT  .  EQ.  2)  THEN 

CALL  GRTRK(  AINRAY , APRAY , TARAY , STF ,NUM ,TIMRAY) 
ELSEIF  (INPUT  .  EQ.  3)  THEN 

GOTO  360 
ELSE 

PRINT*, 'INVALID  ENTRY! ' 

GOTO  340 
END  IF 

DISPLAY  DATA 

CALL  DATA( SINC, SA,SE ,TFEA,TFSU,TFMO,TFDRA, SPER, SPI, TDI,TDA,TDE, 
+         TDMM , TDMA , TDLAN , TDH , TDAP , SMM , SMA , SLAN , SH , SAP , SV , SR) 
CALL  ENDPL(O) 

PROMPT  USER  IF  ANOTHER  DISPLAY  TYPE  IS  DESIRED 

PRINT*, 'WOULD  YOU  LIKE  ANOTHER  PLOT  USING  THE  SAME  ORBITAL' 

PRINT*, 'PARAMETERS  AND  DATA: ' 

PRINT*, 'ENTER  "Y"  OR  "N"  :' 

READ*,YORN 

PRINT*, YORN 


ORB23450 
ORB23460 
ORB23470 
0RB23480 
0RB23490 
ORB23500 
ORB23510 
ORB23520 
ORB23530 
ORB23540 
ORB23550 
ORB23560 
ORB23570 
0RB23580 
ORB23590 
0RB23600 
0RB23610 
0RB23620 
0RB23630 
0RB23640 
ORB23650 
0RB23660 
0RB23670 
ORB23680 
ORB23690 
ORB23700 
0RB23710 
ORB23720 
ORB23730 
0RB23740 
0RB23750 
ORB23760 
ORB23770 
0RB23780 
ORB23790 
ORB23800 
ORB23810 
ORB23820 
ORB23830 
ORB23840 
0RB23850 
0RB23860 
ORB23870 
ORB23880 
0RB23890 
ORB23900 
0RB23910 
ORB23920 
ORB23930 
ORB23940 
ORB23950 
ORB23960 
ORB23970 
ORB23980 
0RB23990 
ORB24000 
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IF  (YORN  .  EQ.  'Y')  THEN 

GOTO  340 
END  IF 
360   RETURN" 
END 

SUBROUTINE  PERIF( RARAY , TARAY ,  NUM ,  P I ,  I NC  ,  LP ,  A ,  E ) 

*  THIS  SUBROUTINE  PLOTS  OUT  THE  RESULTS  OF  THE  PROGRAM  USING  THE 
DISPLAY  FEATURE  ON  THE  MAIN  FRAME. 

REFER  TO  DISSPLA  USERS  GUIDE  FOR  EXPLANATION  OF  DISSPLA 

*  SUBROUTINES. 

REAL  INC, LP 

DIMENSION  TARAY(500) ,RARAY(500) ,RIRAY(500) ,RJRAY(500) ,RKRAY(500) 

I  =  1 

SET  SCALE  OF  AXIS 
RSTEP  =  (A*(l+E))  /  3 
CALL  TEK618 
CALL  RESET(3HALL) 

CALL  SCMPLX 

CALL  PHYSOR(l.  25,4.  ) 

CALL  AREA2D(6.  ,6.  ) 

CALL  MESSAGC  PERIFOCAL  COORDINATE  SYSTEM$  '  ,  100  , 1.  0  ,  6.  5) 

CALL  XNAMEC'XV' ,2) 

CALL  YNAME( 'YW' ,2) 

CALL  XAXANG(90.  0) 

CALL  YAXANG(0.  0) 

CALL  INTAXS 

CALL  POLARf  1.  , RSTEP,  3.  ,3.  ) 

CALL  POLY3 

CALL  NOCHEK 

CALL  CURVE ( TARAY , RARAY , NUM , 1 ) 

CALL  COMPLX 

CALL  HEIGHT(. 2) 

CALL  RESETC COMPLEX') 

CALL  RESETC 'HEIGHT' ) 

CALL  ENDGR(O) 

Display  EARTH  PLOT 

CALL  EARTH1(A,E,INC, PI, RSTEP) 

RETURN 
END 


--■-.J- -•--•- -■--'-  _'.-'., 


.-!--■--'-  -'-.'-.'..'- 


SUBROUTINE  EARTH1( A ,E , INC ,PI , RSTEP) 

THIS  SUBROUTINE  PLOTS  A  VIEW  OF  THE  WORLD,  LOOKING  DOWN  THE  'Z' 

AXIS,  PLACED  ON  THE  ORIGIN.   THE  Latitude  IS  FIXED,  BUT  THE 

LONGITUDE  VARIES  WITH  THE  INCLINATION. 

REFER  TO  DISSPLA  USER"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 

SUBROUTINES 


0RB24010 
0RB24020 
ORB24030 
ORB24040 
ORB24050 
ORB24060 
ORB24070 
ORB24080 
0RB24090 
ORB24100 
ORB24110 
ORB24120 
0RB24130 
0RB24140 
ORB24150 
ORB24160 
ORB24170 
ORB24180 
ORB24190 
ORB24200 
0RB24210 
ORB24220 
0RB24230 
ORB24240 
ORB24250 
ORB24260 
ORB24270 
ORB24280 
ORB24290 
ORB24300 
0RB24310 
ORB24320 
ORB24330 
ORB24340 
ORB24350 
ORB24360 
ORB24370 
0RB24380 
0RB24390 
ORB24400 
0RB24410 
0RB24420 
ORB24430 
0RB24440 
ORB24450 
ORB24460 
ORB24470 
0RB24480 
0RB24490 
ORB24500 
ORB24510 
ORB24520 
ORB24530 
0RB24540 
ORB24550 
ORB24560 
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REAL  INC 

COMMON  IWORK(3800) 

DATA  IWDIM/3800/ 

RE  =  6378. 145 

SCALE  THE  EARTH  PLOT  AND  CENTER  ON  THE  ORIGIN 

SCFAC  =  RE/RSTEP 

SCFAC2  =  SCFAC  *  2.  0 

XPHS  =  1.  25  +  3.  0  -   SCFAC 

YPHS  =  4. 0  +  3. 0  -   SCFAC 

YPOLE  =  90  -  (INC  *  180  /  PI) 

IF(YPOLE  .GT.  90)  THEN 

YPOLE  =  YPOLE  -  90 
END  IF 

YORIG  =  YPOLE  -  90 
YMAX  =  YPOLE  +90 
CALL  RESET(3HALL) 
CALL  PHYSOR(XPHS.YPHS) 
CALL  PRO JCT( ' LAMBERT  EQ/AREA' ) 
CALL  MAPOLE( 0.0, YPOLE) 
CALL  AREA2D  (SCFAC2,SCFAC2) 
CALL  THKFRM  (0.02) 

CALL  GRAF(-90. ,30. ,90. , YORIG, 30. ,YMAX) 
CALL  FRAME 

CALL  MAPFIL('MAPDTA') 
CALL  LBLANK('LAND' ,IWDIM) 
CALL  GRID(1,1) 
CALL  LB LANK( 'WATER' ,IWDIM) 
CALL  DASH 
CALL  GRID(1,1) 
CALL  RESET  ( 'DASH' ) 
CALL  ENDGR(O) 
RETURN- 
END 

SUBROUTINE  GRTRK( AINRAY,APRAY,TARAY,TF,NUM,TIMRAY) 

DIMENSION  AINRAY(500) ,APRAY(500) ,TARAY(500) , 
+   ELARAY(500),ELORAY(500),TLONG(500),TLAT(500),TIMRAY(500) 

RE  =  6. 3782E+03 
EROT  =  7. 292115856E-05 
STF  =  (TF) 
I  =  1 

*     LOAD  ARRAYS  WITH  LATITUDE  AND  LONGITUDE 
410   IF  (I  . LE.  NUM)  THEN 

X  =  RE*COS(APRAY(I))*COS(TARAY(I))-R£*SIN(APRAY(I))* 
+        SIN(TARAY(I)) 

Y  =  RE*COS(AINRAY(I))*SIN(APRAY(I))*COS(TARAY(I))  + 
+       RE*COS(AINRAY( I ) ) *COS( APRAY( I) )*SIN(TARAY( I)) 

Z  =  RE*SIN(AINRAY(I))*SIN(APRAY(I))*COS(TARAY(I))  + 
+       RE*SIN(AINRAY(  I)  )*COS(  APRAY(  I) )*SIN(TARAY(  I) ) 


ORB24570 
ORB24580 
ORB24590 
0RB24600 
0RB24610 
ORB24620 
ORB24630 
ORB24640 
ORB24650 
ORB24660 
ORB24670 
ORB24680 
ORB24690 
ORB24700 
ORB24710 
ORB24720 
ORB24730 
ORB24740 
ORB24750 
ORB24760 
0RB24770 
ORB24780 
ORB24790 
ORB24800 
ORB24810 
ORB24820 
0RB24830 
ORB24840 
ORB24850 
0RB24860 
ORB24870 
ORB24880 
ORB24890 
ORB24900 
0RB24910 
0RB24920 
ORB24930 
ORB24940 
ORB24950 
ORB24960 
ORB24970 
ORB24980 
ORB24990 
ORB25000 
ORB25010 
ORB25020 
ORB25030 
ORB25040 
ORB25050 
ORB25060 
ORB25070 
ORB25080 
ORB25090 
0RB25100 
ORB25110 
ORB25120 
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V- 


CALCULATE  LATITUDE 

ELARAY(I)  =  (ASIN(Z/RE))  *  (180/3.14159) 

TRAP  'X'  AND  'Y'  FOR  ARCTAN  IN  CALCULATING  LONGITUDE 
IF((Y  .  LE.  10)  .AND.  (Y  .  GE.  0.0))  THEN 

Y  =  10. 

ELSEIF  ((Y  .  GE.  -10).  AND.  (Y  .  LE.  0.0))  THEN 

Y  =  -10. 
END  IF 

IF((X  .LE.  10)  .AND.  (X  .  GE.  0.0))  THEN 

X  =  10. 
ELSEIF  ((X  .GE. -10).  AND.  (X  .  LE.  0.0))  THEN 

X  =  -10. 
END  IF 


*  CALCULATE  LONGITUDE 

ELORAY(I)  =  (ATAN2(Y,X)  -  (EROT*TIMRAY( I ) ) ) 

!  MODIFY  LONGITUDES  TO  (  -180  TO  180) 

420     IF  (ELORAY(I)  . LT.  -180)  THEN 

ELORAY(I)  =  ELORAY(I)  +  360 
GOTO  420 
END  IF 
1  =  1  +  1 
GOTO  410 
END  IF 

*  SET  DISSPLA 
CALL  TEK61S 

CALL  RESET(3HALL) 

CALL  YAXANG  (0.  ) 

CALL  PHYSORCl.  0,6.  0) 

CALL  XNAMEC  '  ,1) 

CALL  YNAMEC  '  ,1) 

CALL  AREA2D(7.  5,3.  75) 

CALL  HEADIN  ('GROUND  TRACKS' , 100, 1. 5, 1) 

CALL  SCMPLX 

CALL  MAPGRC-180.  ,90.  ,180.  ,-90.  ,30.  ,90.  ) 

CALL  GRID  (1,1) 

CALL  MAPFIL  ( 'MAPDTA' ) 

I  =  1 

*  IGNORE  Boundary  POINTS 

430   IF  ((ELORAY(I)  .  LT.  -175)  .OR. 
+    (ELORAY(I)  .GT.  175)  .OR. 
+    (ELARAY(I)  .  LT.  -85)  .OR. 
+    (ELARAY(I)  .GT.  85))  THEN 
1  =  1  +  1 
GOTO  430 
END  IF 

ITEMP  =  1 

'<  LOAD  FIRST  POINT  OF  NEW  PLOT  SEGMENT 

IF  (I  .  LE.  NUM)  THEN 


(180/3. 14159) 


ORB25130 
0RB25140 
ORB25150 
ORB25160 
ORB25170 
ORB25180 
ORB25190 
ORB25200 
ORB25210 
ORB25220 
ORB25230 
ORB25240 
ORB25250 
ORB25260 
ORB25270 
ORB25280 
ORB25290 
ORB25300 
ORB25310 
ORB25320 
ORB25330 
ORB25340 
ORB25350 
ORB25360 
ORB25370 
ORB25380 
ORB25390 
ORB25400 
ORB25410 
ORB25420 
ORB25430 
0RB25440 
ORB25450 
ORB25460 
ORB25470 
0RB25480 
ORB25490 
ORB25500 
ORB25510 
ORB25520 
ORB25530 
ORB25540 
ORB25550 
ORB25560 
ORB25570 
ORB25580 
ORB25590 
ORB25600 
ORB25610 
ORB25620 
ORB25630 
ORB25640 
ORB25650 
ORB25660 
ORB25670 
ORB25680 
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TLONG( ITEMP)  =  ELORAY(I) 
TLAT(ITEMP)  =  ELARAY(I) 
1  =  1  +  1 
IF  (  I  . GE.  NUM)  THEN 

CALL  POLYS 

CALL  CURVE(TLONG,TLAT,ITEMP,l) 
END  IF 
END  IF 

LOAD  SECOND  POINT  IN  LINE  SEGMENT 
IF  (I  . LE.  NUM)  THEN 
ITEMP  =  ITEMP  +  1 
TLONG(ITEMP)  =  ELORAY(I) 
TLAT( ITEMP)  =  ELARAY(I) 
1  =  1  +  1 

IF  (  I  .GE.  NUM)  THEN 
CALL  P0LY3 
CALL  NOCHEK 

CALL  CURVE (TLONG,TLAT, ITEMP, 1) 
END  IF 
END  IF 

LOOP  UNTIL  SEGMENT  REACHES  EDGE  OR  NO  MORE  POINTS 

IF  (I  .  LE.  NUM)  THEN 

BOTH  LAT  AND  LONG  INCREASING 
IF((ELORAY(I  -  2)  . LE.  ELORAY( I  -  1))  .AND. 
+      (ELARAYCI  -  2)  .  LE.  ELARAY( I  -  1)))  THEN 

IF((ELORAY(I)  .LT.  -170)  .OR. 
+  (ELARAYCI)  .LT.  -80))  THEN 

CALL  POLY3 
CALL  NOCHEK 

CALL  CURVE (TLONG,TLAT, ITEMP, 1) 
GOTO  430 
ELSE 

ITEMP  =  ITEMP  +  1 
TLONG( ITEMP)  =  ELORAY(I) 
TLAT( ITEMP)  =  ELARAY(I) 
END  IF 

BOTH  LAT  AND  LONG  DECREASING 

ELSEIF((ELORAY(I  -  2)  .  GT.  ELORAY( I  -  1))  .AND. 
+  (ELARAYCI  -  2)  . GT.  ELARAY( I  -  1)))  THEN 

IF((ELORAY(I)  .GT.   170)  .OR. 
+  (ELARAY(I)  .GT.   80))  THEN 

CALL  POLY3 
CALL  NOCHEK 

CALL  CURVE ( TLONG , TLAT , ITEMP , 1 ) 
GOTO  430 
ELSE 

ITEMP  =  ITEMP  +  1 
TLONG( ITEMP)  =  ELORAY(I) 
TLATC ITEMP)  =  ELARAY(I) 
END  IF 

LAT  INCREASING,  LONG.  DECREASING 


ORB25690 
ORB25700 
ORB25710 
ORB25720 
ORB25730 
ORB25740 
ORB25750 
ORB25760 
ORB25770 
ORB25780 
ORB25790 
ORB25800 
ORB25810 
ORB25820 
ORB25830 
ORB25840 
ORB25850 
ORB25860 
ORB25870 
ORB25880 
ORB25890 
ORB25900 
ORB25910 
ORB25920 
ORB25930 
ORB25940 
ORB25950 
ORB25960 
ORB25970 
ORB25980 
ORB25990 
ORB26000 
ORB26010 
ORB26020 
ORB26030 
ORB26040 
ORB26050 
ORB26060 
ORB26070 
ORB26080 
ORB26090 
ORB26100 
ORB26110 
ORB26120 
ORB26130 
ORB26140 
ORB26150 
ORB26160 
ORB26170 
ORB26180 
ORB26190 
ORB26200 
ORB26210 
ORB26220 
ORB26230 
0RB26240 
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+ 
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ELSEIF((ELORAY(I    -    2)    . 

(ELARAYd    -   2)    . 

IF((ELORAY(I)    .  GT. 

( ELARAYd)    •  LT. 

CALL   P0LY3 

CALL 

CALL 

GOTO 


GT.    ELORAY(I 

LE.    ELARAY(I 

170)    .OR. 

-SO))    THEN 


-  1))    .AND. 

-  1)))    THEN 


NOCHEK 
CURVE (TLONG, 
430 


TLAT,ITEMP,1) 


ELSE 


ITEMP  =  ITEMP  +  1 
TLONG(ITEMP)  =  ELORAY(I) 
TLAT(ITEMP)  =  ELARAY(I) 


+ 


ENDIF 

LAT.  DECREASING,  LONG. 

ELSEIF((ELORAY(I    -   2)    . 

(ELARAYd    -   2)    . 

IF((ELORAY(I)    .LT. 

(ELARAYd)    -GT. 

CALL  P0LY3 

CALL 

CALL 

GOTO 


INCREASING 
LE.    ELORAYd 
GT.    ELARAYd 
-170)    .OR. 
80))   THEN 


-  1))  .AND. 

-  1)))  THEN 


NOCHEK 

CURVE ( TLONG, TLAT, ITEMP, 1) 

430 


ELSE 


ITEMP  =  ITEMP  +  1 
TLONG(ITEMP)  =  ELORAYd ) 
TLAT(ITEMP)  =  ELARAYd ) 


END  1 1 

ENDIF 

IF(  I  .  EQ. 

NUM)  THEN 

CALL 

POLY3 

CALL 

NOCHEK 

CALL 

CURVE (TLONG, 

,TLAT, 

ITEMP, 

,1) 

ENDIF 

1  =  1  +  1 

GOTO  440 
ENDIF 

CALL  POLY3 

CALL  NOCHEK 

CALL  CURVE (TLONG, TLAT, ITEMP, 

,1) 

CALL  COMPLX 

CALL  HEIGHT(.  2) 

CALL  THKFRM  (0. 03) 

CALL  FRAME 

CALL  RESETC  COMPLX') 

CALL  RESET( 'HEIGHT') 

CALL  ENDGR  (0) 

RETURN 

END 


ntftihc* 


tffwifyntitrnirfKfeTt 


:-.\--.V 


ORB26250 
ORB26260 
ORB26270 
ORB26280 
ORB26290 
ORB26300 
ORB26310 
ORB26320 
ORB26330 
ORB26340 
ORB26350 
ORB26360 
ORB26370 
ORB26380 
ORB26390 
ORB26400 
ORB26410 
ORB26420 
0RB26430 
ORB26440 
ORB26450 
ORB26460 
ORB26470 
ORB26480 
0RB26490 
ORB26500 
ORB26510 
ORB26520 
ORB26530 
0RB26540 
ORB26550 
ORB26560 
ORB26570 
ORB26580 
ORB26590 
ORB26600 
ORB26610 
ORB26620 
ORB26630 
ORB26640 
ORB26650 
ORB26660 
0RB26670 
ORB26680 
ORB26690 
ORB26700 
ORB26710 
ORB26720 
ORB26730 
ORB26740 
ORB26750 
ORB26760 
ORB26770 
ORB26780 
ORB26790 
ORB26800 
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SUBROUTINE  DATA( I ,A,E ,TFEA,TFSU,TFMO,TFDRA,PER,PI ,TDI ,TDA,TDE , 
+  TDMM , TDMA , TDLAN , TDH , TDAP , MM , MA , LAN , H , AP , V , R ) 

THIS  SUBROUTINE  Displays  THE  ORBITAL  DATA  FOR  BOTH  THE  PERIFOCAL 
AND  THE  GROUND  TRACK  PLOTS. 

REFER  TO  DISSPLA  USER"S  MANUAL  FOR  EXPLANATION  OF  DISSPLA 
SUBROUTINES 

REAL  I, MM, MA, LAN 

MU  =  3. 986012E+05 

CALCULATE  THE  AVERAGE  FORCES  FROM  THE  TOTAL  MAGNITUDE  OF 

*  FORCE  CHANGES 
AVGFE  =  TFEA/50.  0 
AVGFS  =  TFSU  /  50.  0 
AVGFM  =  TFMO  /  50.  0 
AVGFD  =  TFDRA  /  50.  0 

*  CALCULATE  ORBITAL  ELEMENTS  IN  Usable  UNITS 
PERH  =  PER/3600 

DI  =  I  *  (180.0/PI) 
DLAN  =  LAN  *  (180.0/PI) 
DAP  =  AP  *  (180.0/PI) 

CALCULATE  Average  CHANGE  IN  ELEMENTS  FOR  ONE  PERIOD 

AVGDI  =  TDI  /  50.  0 

AVGDA  =  TDA  /  50.  0 

AVGDE  =  TDE  /  50.  0 

AVGDMM  =  TDMM  /  50.  0 

AVGDMA  =  TDMA  /  50.  0 

AVGLAN  =  TDLAN  /  50.  0 

AVGDH  =  TDH  /  50.0 

AVGDAP  =  TDAP  /  50. 0 

CALCULATE  RADIUS'S  AND  VELOCITIES 

ENR  =  ((V**2)/2)  -  (MU/R) 

RP  =  A*(l  -  E) 

RA  =  A*(l  +  E) 

VP  =  SQRT(2*(ENR  +  (MU/RP))) 

VA  =  SQRT(2*(ENR  +  (MU/RA))) 


SET  DISSPLA 

CALL  RESET(3HALL) 

CALL  SCMPLX 

CALL  PHYSOR(0.  0,0.  0) 

CALL  AREA2D(8.  5,4.  0) 

PRINT  DATA 

CALL  MESSAGC'I  =  S1  ,  100,0. 25 ,3. 67) 

CALL  REALNO(DI,3,  ABUT' ,'ABUT') 

CALL  MESSAG( '  DEG. $ ' , 100 , ' ABUT' , ' ABUT' ) 

CALL  MESSAGC   A  =  $',  100  ,'  ABUT"  ,'  ABUT'  ) 

CALL  REALNO( A , 1 , ' ABUT ' , ' ABUT ' ) 

CALL  MESSAG('  KMS ', 100 ,' ABUT' ,' ABUT' ) 


ORB26810 
ORB26820 
ORB26830 
ORB26840 
ORB26850 
ORB26860 
ORB26870 
ORB26880 
ORB26890 
ORB26900 
ORB26910 
ORB26920 
ORB26930 
ORB26940 
ORB26950 
ORB26960 
ORB26970 
ORB26980 
ORB26990 
ORB27000 
ORB27010 
ORB27020 
ORB27030 
ORB27040 
ORB27050 
ORB27060 
ORB27070 
ORB27080 
ORB27090 
ORB27100 
ORB27110 
ORB27120 
ORB27130 
ORB27140 
ORB27150 
ORB27160 
ORB27170 
ORB27180 
ORB27190 
ORB27200 
ORB27210 
ORB27220 
ORB27230 
ORB27240 
ORB27250 
ORB27260 
ORB27270 
ORB27280 
ORB27290 
ORB27300 
ORB27310 
ORB27320 
ORB27330 
ORB27340 
ORB27350 
ORB27360 
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CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

MESSAGC 

+ 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

R£ALNO( 

CALL 

MESSAG( 

CALL 

REALNOC 

CALL 

MESSAG( 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

+$' ,100,1. 0,1 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

REALNOC 

CALL 

MESSAGC 

CALL 

MESSAGC 

CALL 

REALNOC 

'  E  =  S' ,100, 'ABUT' ,'ABUT') 

E, 3, 'ABUT' , 'ABUT' ) 

'   PER  =$' ,100, 'ABUT' ,'ABUT1) 

PERH,2,'ABUT' ,'ABUT* ) 

'  HOURS$' ,100, 'ABUT' ,' ABUT' ) 

'AVERAGE  RATE  OF  CHANGE  OF  ELEMENTS  PER  SECOND  $', 
00,1.  0,3.  0) 

'DI/DT  =  $' ,100,0. 25,2.67) 

AVGDI,-2,'ABUT' ,'ABUT') 

'   DA/DT  =  $' ,100 ,'ABUT' ,'ABUT1) 

AVGDA,-2,'ABUT' ,'ABUT' ) 

1   DE/DT  =$' ,100, 'ABUT' ,'ABUT') 

AVGDE, -2 ,'ABUT' ,'ABUT' ) 

'DMM/DT  =  $' ,100,0. 25,2. 33) 

AVGDMM,-2,  ABUT' ,'ABUT') 

'   DMA/DT  =  $' ,100, 'ABUT' ,'ABUT' ) 

AVGDMA,-2,'ABUT' ,'ABUT') 

'   DLAN/DT  =  $' ,100, 'ABUT' ,'ABUT' ) 

AVGLAN , - 2 , ' ABUT ' , ' ABUT ' ) 

'DH/DT  =  $' ,100,0. 25,2. 00) 
AVGDH,-2,'ABUT' ,'ABUT' ) 
'   DAP/DT  =$' ,100, 'ABUT' ,'ABUT') 
AVGDMA , - 2 , ' ABUT ' , ' ABUT ' ) 

'AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S**2) 

.67) 

'EARTH  =  $'  ,100,0.  10,1. 33) 

AVGFE.-l, 'ABUT' ,'ABUT' ) 

'  MOON  =  $' ,100, 'ABUT' ,'ABUT') 

AVGFM,-1,'ABUT' ,'ABUT' ) 

'  SUN  =  $' ,100, 'ABUT' ,'ABUT') 

AVGFS,-1,'ABUT' ,'ABUT1 ) 

'  DRAG  =  $' ,100, 'ABUT' ,'ABUT') 

AVGFD,-1,'ABUT' ,'ABUT1) 

'PERIGEE$'  ,100,2. 75,1. 0) 

Apogee$ ', 100 ,' ABUT' ,' ABUT' ) 

'RADIUS  (KM)S'  ,100,0.25,0.67) 
'RP  =$' ,100,2. 75,0. 67) 
RP,1,'ABUT' ,'ABUT' ) 

$' ,100, 'ABUT' ,'ABUT') 

'   RA  =$' ,100, 'ABUT' ,'ABUT' ) 
RA,1,'ABUT' ,'ABUT') 

'VELOCITY  ( KM/SEC) $' ,100,0.  25,0.  33) 
*VP  =$'  ,100,2. 75,0. 33) 
VP,2,'ABUT' ,'ABUT' ) 

$'  ,100, 'ABIT' ,'ABUT' ) 
'   VA  =$' ,100, 'ABUT' ,'ABUT' ) 
VA,2, 'ABUT' ,'ABUT' ) 


ORB27370 
ORB27380 
ORB27390 
ORB27400 
ORB27410 
ORB27420 
ORB27430 
0RB27440 
ORB27450 
ORB27460 
ORB27470 
ORB274S0 
ORB27490 
ORB27500 
ORB27510 
ORB27520 
ORB27530 
ORB27540 
ORB27550 
ORB27560 
ORB27570 
ORB27580 
ORB27590 
ORB27600 
ORB27610 
ORB27620 
ORB27630 
ORB27640 
ORB27650 
ORB27660 
ORB27670 
ORB27680 
ORB27690 
ORB27700 
ORB27710 
ORB27720 
ORB27730 
ORB27740 
ORB27750 
ORB27760 
ORB27770 
ORB27780 
ORB27790 
ORB27800 
ORB27810 
ORB27820 
ORB27830 
ORB27840 
ORB27850 
ORB27860 
ORB27870 
ORB27880 
ORB27890 
ORB27900 
ORB27910 
ORB27920 
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CALL  RESET('COMPLX')  0RB2' 

CALL  ENDGR(O)  ORB2  7 

RETURN  ORB 2 7 

END  ORB 2 7 

ORB  2  7 
ORB  2  7 


APPENDIX  B.     COORDINATE  SYSTEMS 
A.     'UK'  :  GEOCENTRIC  -  EQUATORIAL 


*K 


vernal  equinox 
direclion 


Figure  3.       Geocentric-equatorial  coordinate  system 

The  geocentric-equatorial  system  as  seen  in  Figure  3  has  its  origin  at  the  earth's 
center.  The  fundamental  plane  is  in  the  equator  and  the  positive  X-axis  points  in  the 
vernal  equinox  direction.  The  Z-axis  points  in  the  direction  of  the  north  pole.  This 
system  is  not  fixed  to  the  earth  and  turning  with  it;  rather,  the  geocentric-equatorial 
frame  is  nonrotating  with  respect  to  the  stars  (except  for  precession  of  the  equinoxes) 
and  the  earth  turns  relative  to  it.  Unit  vectors.  I.J,  and  A'  shown  in  Figure  3.  lie  along 
the  X,  Y.  and  Z  respectively.   [Ref.  1:  p. 55] 
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B.     'PQW:  PERIFOCAL 


Figure  -4.        Perifocal  coordinate  system 

The  perifocal  coordinate  system  has  its  fundamental  plane  in  the  plane  of  the  satel- 
lite's orbit  as  seen  in  Figure  4.  The  coordinate  axis  are  named,  X„,  Yu  and  Z„  .  The 
X„  axis  points  toward  the  perigee;  the  YM  axis  is  rotated  90  degrees  in  the  direction  of 
orbital  motion  and  lies  in  the  orbital  plane;  the  Zu  axis  along  h  completes  the  right- 
handed  perifocal  system.  Unit  vectors  in  the  direction  of  X„,  Y^  and  Z^  are  called  P,  Q 
and  IV  respectively.   [Ref.  1:  p. 57] 
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C.     'RSVV  :  ORBITAL 


z       iS 
W                 '      \      J 

X 

\     r   /   " — -^ 

\    /       ^^^R 

i 

Figure  5. 


Orbital  coordinate  svstem 


(Figure  9.4-1.  Ref.  1) 
The  orbital  coordinate  system  has  its  principle  axis,  R  (unit  vector  r),  along  the  in- 
stantaneous radius  vector,  r  as  seen  in  Figure  5.  The  axis  S  is  rotated  90  degrees  from 
R  in  the  direction  of  increasing  true  anomaly.  The  third  axis,  W,  is  perpendicular  to 
both  R  and  S.  Note  that  this  coordinate  system  is  simply  rotated  v„  from  the  PQW 
perifocal  system.   [Ref.  1:  p.39S] 
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D.     COORDINATE  TRANSFORMATIONS 

The  coordinate  transformations,  for  the  previous  coordinate  systems,  use  angular 
rotations  about  the  axis  to  evaluate  the  transformation  matrix.  The  matrix  elements  ru 
are  calculated,  then  applied  to  the  old  vector  to  get  the  vector  in  the  new  coordinate 
system.   The  following  orbital  elements  are  used: 

Q  =  longitude  of  ascending  node 

co  =  argument  of  perigee 

i  =  inclination 

u0  =  argument  of  latitude 

v0  =  true  anomaly 
The  coordinate  transformations  follow  [Ref.  1:  p.74-83] 

1.  PQWto  UK 

>-,,  =  cos  Q.  cos  co  —  sin  Q  sin  co  cos  i 

rn  =  —  cos  il  sin  co  —  sin  Q.  cos  co  cos  / 

fi3  =  sin  Q.  cos  co 

r2.  =  sin  Q  cos  co  +  cos  Q  sin  co  cos  i 

r22  =  —  sin  Q  sin  co  +  cos  Q  cos  co  cos  / 

r:3  =  —  cos  Q  sin  i 

/•;,  =  sin  co  sin  i 

ri:  —  cos  co  sin  / 

;• ..  =  cos  i 

7  =  rnP  + rl2Q  + ralV 
J  =  ,-2lP  +  r220  +  rrJV 
k=r3lP  +  r32Q  +  r33lV 

2.  UK  to  PQW  (inverse  of-1) 

P  =  rj  +  r2]J  +  r3xK 
Q  =  rl27  +  r22J  +  r32K 

IV  =  >\J  +  rnJ  +  rVjK 

3.  UK  to  RSYV 

/-,,  =  cos  Q.  cos  u0  —  sin  Q.  sin  u„  cos  i 
r,2  =  sin  Q  cos  u0  +  sin  u0  cos  Q.  cos  / 
rl3  =  sin  /  sin  wn 

r2,  =  —  cos  LI  sin  w0  —  sin  Q.  cos  u0  cos  i 
r22  =  —  sm  Q  sin  u0  +  cos  Q.  cos  w0  cos  / 
y:3  =  cos  w0  sin  i 

«•_.  =  <;m  O  sin  /' 


'31 

''3: 

=  —  cos 

Mil    1 

Q  sin  / 

r33 
R 

=  COS  / 

r,.J  +  r4 

ik 

S 

-r3I7  + 

rj  +  > 

}k 

■J 

SO 


4.  RSW  to  UK  (inverse  of -3) 

7  =  ruR  +  r3lS  +  /s,M' 
J  =  r12^  +  r22S  +  rjjjf 
A'  =  r13/?  -  r2iS  +  rH  H ' 

5.  PQW  to  RSW 

r„  =  cos  v„ 
ru  =  sin  v0 
r13  =  0.0 

/•:,  =  COS  v„ 
r,3  =  0.0 
/•31  =  0.0 
r32  =  0.0 

>33=  1.0 

tf  =  rnP  +  /-12£  +  r13J?' 
5  =  r:iP  +  r2:0  +  r2ilV 
IV  =  hlP  +  rnQ  +  rJV 

6.  RSW  to  PQW  (inverse  of  #5) 

P  =  >uR-rr2lS  +  r:JV 
0  =  rnR  +  r^S  +  r22W 
W  =  r13R  +  raS  +  raW 
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APPENDIX  C.     ORBITAL  ELEMENTS 

The  user  is  assumed  to  be  studying  orbital  mechanics  and  should  understand  the 
orbital  elements  and  how  to  calculate  them.  A  brief  description  of  the  elements  and  the 
equations  used  to  calculate  the  elements  follow.  For  a  detailed  explanation  of  the  ele- 
ments and  the  equations  to  calculate  them  refer  to  Chapters  1  and  2  of  reference  1. 
Figure  6  on  page  83  shows  the  orbital  elements  in  the  Geocentric-Equatorial  and 
perifocal  coordinate  system. 

1.   Angular  Momentum  (h): 

The  specific  angular  momentum  is  a  constant  of  the  motion  of  the  satellite. 

defined  as  h  =  rxv. 

h  =  7x7  =  htI  +  hjJ  +  hkK 

k  =  >jvk  -  rkvj 

hi  =  r,.v.  —  r[Vk 

ll:  =   J';'.";  —   /•;V: 


h  =  N    hT  +  hj  +  hi 

2.   Node  Vector  (n): 

The  node  vector  is  a  vector  pointing  along  the  line  of  nodes  in  the  direction 

of  the  ascending  node. 

7  =  Kxh  =  -II  I  +  hj 


n  =  N  hj  ~  lh 

3.   Semi-latus  rectum  (p): 

The  semi-latus  rectum  is  a  geometric  constant  of  the  conic  section. 

/;- 
P  =  — 

A.   Eccentricity  (e): 

The  eccentricity  is  a  constant  defining  the  shape  of  the  conic  orbit. 

e  =  —  [(v2-  —  )'>:-{77)7] 


5.   Semi-major  axis  in): 


S2 


fr-rilt    «•.«/*„  . 


Figure  6.      Orbital  elements 

The  semi-major  axis  is  a  constant  defining  the  size  of  the  orbit. 

(1  -e2) 

a  = - 

6.   Inclination  (i): 

The  inclination  is  the  angle  between  the  K'  unit  vector  in  the  TJK'  system  and 
the  angular  momentum  vector,  'h'. 


S3 


-l    h  K  -l  -  hk 

i  =  COS     (  — ; —  )  =  COS     (  — — 
1     h  n 


Longitude  of  ascending  node  {Q.Y. 

The  longitude  of  the  ascending  node  is  the  angle  in  the  fundamental  plane  , 
between  the  T  unit  vector  and  the  point  where  the  satellite  crosses  through  the 
fundamental  plane  in  a  northerly  direction  (ascending  node)  measured  counter- 
clockwise when  viewed  from  the  north  side  of  the  fundamental  plane. 


1  1: 

n  =  cos-'(  — ) 


8.  Argument  of  perigee  (co): 

The  argument  of  perigee  is  the  angle  in  the  plane  of  the  satellite's  orbit,  be- 
tween the  ascending  node  and  the  perigee  point,  measured  in  the  direction  of  the 
satellite's  motion. 


-1/  in 


')  =  cos     t— '  =  co 


(nft  +  nje,) 


HC 


9.   True  anomaly  at  epoch  (v0): 

The  true  anomaly  at  epoch  is  the  angle  in  the  plane  of  the  satellite's  orbit,  be- 
tween perigee  and  the  position  of  the  satellite  at  a  particular  time.  t0,  called  the 
"epoch". 


-1/  e  i 

Vn  =  COS      ( 


er    ' 

10.  Argument  of  latitude  (u„): 

The  argument  of  latitude  is  the  angle  in  the  plane  of  the  orbit,  between  the 
ascending  node  and  the  radius  vector  to  the  satellite  at  time  iv. 

-1/  nr   x 
«n  =  cos    ( 


nr 


11.  Longitude  of  perigee  (IT): 

The  longitude  of  perigee  is  the  angle  from  T  to  perigee  measured  eastward  to 
the  ascending  node  and  then  in  the  orbital  plane  to  perigee. 

n  =  a  +  co 

12.  True  longitude  at  epoch  (/„): 

The  true  longitude  at  epoch  is  the  angle  between  T  and  r0  (the  radius  vector 
to  the  satellite  at  i0  measured  eastward  to  the  ascending  node  and  then  in  the  orbital 
plane  to  r0. 


L  =  co  -I-  Q.  +  v 


0 


13.  Period  (per): 

The  period  is  the  time  the  for  the  satellite  to  complete  one  orbit. 


Per  =  2  -^— 


14.  Eccentric  anomalv  (LA): 
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The  eccentric  anomaly  is  the  angle  between  the  perigee  and  a  position  on  an 
auxiliary  circle  circumscribed  about  the  ellipse  where  a  perpendicular  line  to  the 
major  axis  has  been  extended  from  the  epoch  location  of  the  satellite  to  the  auxil- 
iary circle. 

_,     e  +  cos(v) 

EA  =  cos 


1  -f  e  cos(v) 


15.  Mean  motion  (n'): 

The  mean  motion  is  defined  below: 


/  u 


V 

16.  Mean  anomaly  (MA): 

The  mean  anomaly  is  defined  below: 

MA  =  n'ii-  T)  =  EA  -  e  smiEA) 

17.  Time  of  flight  (TF): 

The  time  of  flight  is  the  elapsed  time  from  when  the  satellite  was  at  perigee  to 
the  current  epoch. 


3 


I  -  T)  =  yj  —(EA-e  s'm(EA)) 
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APPENDIX  D.     SAMPLE  ORBITS 

To  demonstrate  the  capabilities  of  the  program,  a  variety  of  orbital  plots  will  follow: 
1.    Low  earth  orbit  (LEO). 


PERIFOCAL  COORDINATE  SYSTEM 


I  =    45.000  DEC.    A 


7*44.4 

7222.2  KM  E 


0.100    PER  «    1.70  HOURS 


AVERAGE  RATE  OP  CHANGE  OP  ELEMENTS  PER  SECOND 
DI/DT  =    0.00    DA/DT-    0.00    DE/DT  «    0.00 
DMM/DT«    0.00    DMA/DT-    0.00    DLAN-DT »    0.00 
DH/DT=    0.00    DAP/DT-    0.00 

AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S"2> 
EARTH  =    0.0  MOON  ■    0.0  SUN  *    0.0  DRAG  -    0.0 

PERIGEE  Apogee 

RADIUS  (KM)  RP  -  6500.0  RA  -  7044.4 

VELOCITY  (KM/SEC)      VP  «  8.21  VA  «  6.72 


Figure  7.      Unperturbed  Low  Earth  Orbit  (LEO) 


Figure  7  shows  the  perifocal  plot  of  a  satellite  in  an  unperturbed  low  earth 
orbit  (LEO).   The  initial  parameters  of  the  orbit  were  entered  as  follows: 
radius  of  perigee  (RP)  =  6500  km 
eccentricity  (e)  =  0.1 
inclination  (i)  =  45  degrees. 
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PERIFOCAL  COORDINATE  SYSTEM 


T»07- 

I  =    44.998  DEC.    A  -    7203.3  KM  E  =    0098    PER 


1.69  HOURS 


AVERAGE  RATE  OF  CHANGE  OF  ELEMENTS  PER  SECOND 
DI/DT=    4.2OM0"*    DA/DT  =    e^S'lO"*    DE/DT  *    1.00'10~* 
DMM/DT=    1.80#10"'    DMA/DT-    9.21 '10"*    DLAN/DT  =    9.48M0"* 
DH/DT  =    5.83'10*"    DAP/DT=    9-21*10"* 

AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S~ 2) 
EARTH  -    9.6M0**  MOON  «=    9.4'10"**  SUN  -    4-3"10~**  DRAG  *    1.4*10" 


RADIUS  (KM) 
VELOCITY  (KM/SEC) 


PERIGEE 
RP  -  6499.6 
VP  =  8.20 


Apogee 
RA  «■  7907.0 
VA  «=  6.74 


Figure  8. 


Perturbed  Low  Earth  Orbit  (LEO) 


With  perturbing  forces  applied  to  the  previous  LEO,  the  drag  force  will  be  the 
dominate  perturbing  force.  The  drag  will  act  as  a  negative  velocity  change  applied 
in  the  area  of  perigee,  with  the  result  of  decreasing  the  semi-major  axis  length,  this 
in  effect  will  decrease  the  eccentricity  of  the  orbit,  as  can  be  seen  by  comparing  the 
orbital  data  of  the  unperturbed  LEO  in  Figure  7  on  page  86  with  the  orbital  data 
of  the  perturbed  LEO  in  Figure  S. 
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2.   Circular  orbit. 


GROUND  TRACK 


•C*  H 


I  =    59.967  DEC.    A  *    6986.1  KM  E  -    0000    PER  =    1.61  HOURS 

AVERAGE  RATE  OF  CHANGE  OF  ELEMENTS  PER  SECOND 
D1/DT  =    4.02*10"*    DA/DT-    9.74*10"*    DE/DT «    5.93*10"" 
DMM/DT  =    2-26*10"*    DMA/DT  =    IJX'IO"    DLAN/DT  =    7.32*10"* 
DH/DT-    5.86*10"*    DAP/DT-    1.52*10"" 

AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S*"2) 
EARTH  «    IJZ'10"  MOON  -    8.8*10"  SUN  =    4.3*10""  DRAG  «    3.0'10"" 


PERIGEE 
RADIUS  (KM)  RP  -  6984.4 

VELOCITY  (KM/SEC)      VP  =  7.56 


Apogee 
RA  »  6987.8 
VA  -  7.55 


Figure  9.      Circular  Orbit 

An  example  of  the  plot  of  the  ground  track  of  a  sequence  of  three  60  degree 
inclined  perturbed  circular  orbits  with  a  radius  of  7000  km  is  shown  in  Figure  9. 
The  sequence  of  orbits  displays  the  precession  of  the  orbit  around  the  earth. 
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3.   Transfer  orbit. 


PERIFOCAL  COORDINATE  SYSTEM 


ttOtl  J-, 


XW 


I  =    0.000  DEC.    A  =    14510.8  KM  E  ■    0.518    PER  -    4.83  HOURS 

AVERACE  RATE  OF  CHANCE  OF  ELEMENTS  PER  SECOND 
DI/DT  =    0.00    DA/DT  =    0.00    DE/DT  =    0.00 
DMM/DT  =    0.00    DMA/DT  =    0.00    DLAN/DT  «=    0.00 
DH/DT  «    0.00    DAP/DT  =    0.00 

AVERACE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KU/S"Z) 
EARTH  =    0.0  MOON  -    0.0  SUN  •=    0.0  DRAG  =    0.0 

PERIGEE  Apogee 

RADIUS  (KM)  RP  -  7000.0  RA  «  22021.5 

VELOCITY  (KM/SEC)      Vp  =  fi.30  VA  -  2.95 


Figure   10.      Transfer  Orbit 


The  transfer  orbit  between  a  circular,  equatorial  LEO  and  a  molniya  orbit 
(high  eccentric  orbit)  is  shown  in  Figure  10.  A  velocity  increase  of  1.75  km/s  was 
applied  at  the  perigee  to  simulate  a  perigee  kick  to  boost  the  satellite  into  the 
molniya  orbit.  A  similar  velocity  change  could  then  be  applied  at  apogee  to  create 
a  high  altitude  circular  orbit,  or  a  negative  velocity  change  applied  at  perigee  could 
be  used  to  bnns  the  satellite  back  to  a  LEO. 
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4.   Geosvnchronous  orbit 


GROUND  TRPCK 


1  «    60.000  DEC.    A  -    42234.6  KM  E  =    0.000    PER  -    23.99  HOURS 

AVERAGE  RATE  OF  CHANCE  OF  ELEMENTS  PER  SECOND 
Dl/DT  =    1.04*10"*    DA/DT=    1.26'10**    DE/DT  =    1.08M0"* 
DMM/DT  =    3.26M0"    DMA/DT  =    3.02MO"    DLAN/DT  ■    LSS'IO"* 
DH/DT  =    3.M'10~*    DAP/DT«    3-02'10"' 

AVERAGE  MAGNITUDE  OF  FORCES  PER  UNIT  MASS  (KM/S*"2) 
EARTH  »    9.3"10~*  MOON  «    4-9'10"*  SUN  -    2.6*10"*  DRAG  «    3.1"10"'< 

PERIGEE  Apogee 

RADIUS  (KM)  RP  ■  42233.9  RA  •  42233.3 

VELOCITY  (KM/SEC)      VP  =  3.07  VA  -  3.07 


Figure   11.      Geosynchronous  Orbit 


The  ground  track  of  a  perturbed  geosynchronous  orbit  inclined  60  degrees  is 
shown  in  Figure  11.  The  orbit  displays  the  figure  eight  typical  with  inclined 
geosvnchronous  orbits. 
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